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Abstract 


A  simulation  model  specifically  designed  to  determine  the 
probable  hydrologic  changes  resulting  from  watershed  management 
in  the  Colorado  subalpine  zone  is  described.  The  model  simulates  the 
total  water  balance  on  a  continuous  year-round  basis  and  compiles 
the  results  from  individual  hydrologic  response  units  into  a  "com- 
posite overview"  of  an  entire  drainage  basin.  Preliminary  results  are 
summarized  for  an  8-year  test  period  on  a  667-acre  experimental 
watershed. 
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Hydrologic  Simulation  Model  of 
Colorado  Subalpine  Forest 

Charles  F.  Leaf  and  Glen  E.  Brink 


Leaf  and  Brink  (1973)  have  previously 
described  a  model  for  simulating  snowmelt  in 
central  Colorado  subalpine  watersheds. 
Snowmelt  over  an  area  is  described  in  terms  of 
combinations  of  aspect,  slope,  elevation,  and 
forest  cover  composition  and  density. 

The  hydrologic  model  described  in  this 
report  is  an  expanded  version  of  the  snowmelt 
model.  The  model  has  been  programed  for  the 
CDC  6400  computer  at  Colorado  State  Univer- 
sity. It  is  designed  to  simulate  the  total  water 
balance  on  a  continuous,  year-round  basis,  and 
to  compile  the  results  from  individual 
hydrologic  subunits  into  a  "composite  over- 
view" of  an  entire  watershed.  The  model  has 
been  designed  to  simulate  watershed  manage- 
ment practices  and  their  resultant  effects  on  the 


behavior  of  hydrologic  systems.  The  model 
consists  of  (1)  a  "core"  which  performs  the 
actual  simulation,  and  (2)  peripheral  routines 
which  specify  hydrologic  subunit  parameters, 
obtain  the  input  data,  maintain  continuity 
between  simulation  intervals,  and  output  the 
results. 

Figure  1  schematically  shows  the 
general  flow  of  the  model.  Detailed  flow  chart 
descriptions  of  the  water  balance  routines  and 
pertinent  hydrologic  theory  are  presented  in 
this  report.  Those  routines  which  were  incor- 
porated from  the  snowmelt  model  without 
significant  changes  are  not  discussed  here. 
Complete  descriptions  of  the  unrevised  routines 
are  given  in  Leaf  and  Brink  (1973).  The  routines 
which  were  taken  from  the  snowmelt  model  and 
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Figure  1.  — General  flow  of  hydrologic 
simulation  model. 
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altered  for  compatibility  with  the  water 
balance  model  are  indicated  by  an  asterisk  (*) 
in  the  following  tabulation. 

Discussed  in  this  report 

WATBAL* 

CANVAP 

SNOWVAP 

EVTRAN 

RADBAL* 

SNOWED* 


Discussed  in  Leaf  and  Brink  (1973) 

AFFECTS 
CALIN 
CALOSS 
DIFMOD 

GETREF  (now  PACKREF) 

LINK 

MIXTURE 

RAINED 

Subroutine  WATBAL  (fig.2) 


Evapotranspiration 

Morton  (1971)  points  out  that  "the 
relationship  between  potential  evaporation 
and  regional  (actual)  evaporation  includes  the 
effects  of  hydrologic  and  climatologic  feed- 
back." The  feedback  includes  moisture  supply 
and  the  thermal  and  moisture  characteristics  of 
the  overlying  air,  which  are  influenced  by  the 
actual  evapotranspiration.  This  interaction  in 
turn  has  a  significant  influence  on  the  energy 
available  for  evapotranspiration. 

These  interactions  have  also  been  taken 
into  account  by  Bouchet  (1963),  who  argued 
that  changes  in  regional  and  potential  evapora- 
tion due  to  changes  in  regional  moisture  supply 
are  complementary.  If  the  potential 
evapotranspiration  (ET)  is  computed  from 
regional  climatological  observations  and  utiliz- 
ed in  this  concept,  the  regional  (actual) 
evapotranspiration,  which  is  a  product  of  com- 
plex climatic,  soil  moisture,  and  vegetative 
processes  may  be  estimated. 

Of  the  several  empirical  methods 
available  for  computing  potential 
evapotranspiration,   the  one  developed  by 


Subroutine  AFFECTS  from  the 
snowmelt  model  (Leaf  and  Brink  1973)  was 
expanded  to  include  the  decisions  relating  to 
evapotranspiration,  and  was  renamed  WAT- 
BAL. WATBAL  is  the  primary  routine  in  the 
water  balance  model.  It  receives  input  on  a 
daily  basis,  the  subunit  parameters,  and  all 
state  variables  computed  by  the  peripheral 
routines.  (The  only  links  between  WATBAL 
and  the  peripheral  routines  are  common  block 
/WATRBAL/  and  the  formal  parameters  pass- 
ed at  the  time  of  the  call.) 

Precipitation 

Precipitation  (if  ^ny)  is  classified  as 
discussed  in  AFFECTS  (Leaf  and  Brink  1973), 
and  the  degree  to  which  it  affects  the  energy 
balance  is  calculated. 


Figure  2.  —  Subroutine  WATBAL 
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Hamon  (1961)  appears  to  give  the  best  results 
(Stephens  and  Stewart  1963,  Russell  and 
Boggess  1964,  Takhar  and  Rudge  1970). 
Hamon  formulated  the  expression: 


=  CD^Pj.  [1] 


for    computing    average  potential 
evapotranspiration,  E^^  (inches /day),  for  each 
month  of  the  year 
where 

D  =  possible  sunshine  in  units  of  12  hours, 
=  The  saturated  water  vapor  density  (ab- 
solute  humidity)   at  the  daily  mean 
temperature  in  grams  per  cubic  meter, 
and 

C  =  a  coefficient  (0.0055  according  to  Hamon). 

Hamon's  equation  requires  only  latitude,  con- 
verted to  day  length  (adjusted  for  slope  and 
aspect),  and  mean  temperature,  converted  to 
saturation  vapor  density,  for  computing 
monthly  Ej^. 

Equation  [1]  predicts  the  "average" 
evapotranspiration.  In  the  Colorado  subalpine 
zone,  this  average  is  less  than  half  the  amount 
that  could  occur  under  conditions  of  unlimited 
energy  supply,  assumed  herein  as  potential 
solar  radiation.  Accordingly,  the  coefficient  C 
in  equation  [1]  was  empirically  adjusted  up- 
ward to  obtain  an  expression  for  potential 
evapotranspiration  under  maximum  solar  in- 
put: 

Em  =  C'D^P^.  [2] 

where  C '  is  the  adjusted  coefficient.  The  daily 
potential  evapotranspiration  for  each  of  12 
months  as  derived  by  equation  [2]  is  supplied  as 
a  set  of  parameters  for  each  hydrologic  subunit. 

To  adjust  maximum  daily  evapo- 
transpiration for  available  energy,  the  values 
determined  by  equation  [2]  were  modified  ac- 
cording to  the  expression 


where 

E    =  evapotranspiration  adjusted  for  available 

energy  in  inches/day, 
SW  =  the  observed  daily  shortwave  radiation  in 

langleys, 

P    =  potential  shortwave  radiation  for  the  day 
as  computed  by  Frank  and  Lee  (1966). 


In  the  water  balance  routines,  the  ad- 
justed evapotranspiration  as  derived  above  is 
then  redefined,  depending  on  the  source,  as 
selected  by  the  following  sequence: 

1.  If  snow  is  intercepted  on  the  forest  canopy, 
evaporation  occurs  exclusively  from  that 
source  and  is  computed  by  subroutine  CAN- 
VAP. 

2.  If  the  canopy  is  free  of  snow,  the  next  step  in 
the  source  selection  is  to  determine  if  losses 
result  from  evapotranspiration  (see  sub- 
routine EVTRAN)  or  evaporation  from  the 
snowpack  surface  (see  subroutine  SNOW- 
VAP).  If  evaporation  is  fi:om  the  snow  sur- 
face or  from  intercepted  snow,  control  then 
passes  to  the  radiation  balance  routines.  If  a 
snowpack  exists,  the  radiation  routines 
generate  any  possible  melt  for  input; 
otherwise,  the  only  input  that  can  result  is 
from  a  rain  event.  Subroutine  EVTRAN 
then  calculates  the  evapotranspiration  re- 
quirements, which  are  taken  first  from  the 
input  and,  if  not  satisfied,  from  the  soil 
mantle  storage. 

The  various  methods  of  computing 
evaporation  and  transpiration  are  discussed  in 
the  descriptions  of  the  subroutines  named 
above. 

Once  the  evapotranspiration  re- 
quirements have  been  satisfied,  any  remaining 
input,  either  from  snowmelt  or  rainfall,  is  used 
to  satisfy  the  soil  mantle  recharge  re- 
quirements (see  subroutine  EVTRAN).  When 
field  capacity  is  reached,  the  excess  input  is 
considered  to  be  water  available  for  streamflow 
(generated  runoff. 

As  explained  later  in  this  report,  sub- 
routine RADBAL  includes  a  phase  indicator 
that  determines  which  of  two  methods  is  to  be 
used  to  compute  the  effects  of  the  radiation 
balance.  When  the  seasonal  snowpack  is  com- 
pletely melted,  the  phase  indicator  is  reset  to 
the  "accumulation  phase."  It  remains  at  that 
setting  until  certain  conditions  specified  in 
subroutine  RADBAL  are  met;  it  then  returns  to 
the  "melt  phase"  setting.  Upon  completion  of 
the  water  balance  calculations,  WATBAL 
returns  the  results  and  the  new  values  for  the 
state  variables  to  the  calling  routine.  Here,  they 
are  weighted  according  to  the  percent  of  the 
total  area  occupied  by  the  various  hydrologic 
subunits.  These  weighted  values  are  then 
summed  to  generate  the  watershed  composite. 
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Subroutine  CANVAP  (fig.  3) 

Hoover  and  Leaf  (1967),  Hoover  (1969), 
and  Hoover  (in  press)  ,2  have  discussed  the 
process  and  significance  of  interception  loss  in 
central  Colorado  subalpine  forests.  Field 
studies  indicate  that  mechanical  removal  of 
intercepted  snow  by  wind  is  an  important 
phenomenon.  Accordingly,  wind  effects  were 
considered  in  the  snow  interception  subroutine. 

In  developing  this  portion  of  the  model, 
the  following  assumptions  were  made: 


tion  of  forest  cover  density,  C^,  and 
evapotranspiration  is  adjusted  for  available 
energy  (equation  3)  as  follows: 


where 

Vj,  =  intercepted  snow  evaporation  in  inches. 


The  amount  of  snow  intercepted  varies  ac- 
cording to  forest  cover  type  and  density; 

The  intercepted  snow  rests  on  the  canopy  for 
only  1  day  following  the  day  of  the  snow 
event  because  turbulent  winds  remove  the 
snow  from  the  crowns;  and 

The  residual  intercepted  snow  which  is  not 
vaporized  after  1  day  is  added  to  the  snow- 
pack. 
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The  amount  of  snow  intercepted  by  spruce-fir 
was  assumed  to  vary  as 


NO  SNOW 
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[4] 


dmx 


where 

P-j^f  =  water  equivalent  of  intercepted  snow 
in  inches 

■'■s  =  water  equivalent  input  which  occurs 

as  snow  in  inches, 
Cdmx  =  natural  forest  cover  density,  expressed 

as  a  decimal,  and 
Cjj  =  reduced  forest  cover  density,  as  a 

decimal. 

Interception  in  lodgepole  pine  is  given  by  the 
equation 

C, 

P.    =  0.10       —  [5] 
""■^  dmx 

The  assumed  maximum  amounts  of  snow  in- 
terception are  0.2  and  0.3  inch  for  lodgepole 
pine  and  spruce-fir,  respectively.  Snowfall  in- 
puts which  exceed  the  above  values  are  added 
to  the  snowpack. 

Vaporization  of  intercepted  snow  on 
foliage  surfaces  is  assumed  to  vary  as  a  func- 

'  Hoover,  Marvin  D.  Snow  interception  and  redistribu- 
tion in  the  forest.  Third  Int.  Seminar  for  Hydrol.  Professors 
[Purdue  Univ.,  Lafayette,  Ind.,  July  1971]  (in  press). 


Figure  3.  —  Subroutine  CANVAP. 


If  equation  [6]  yields  a  value  which  is  less  than 
the  water  equivalent  of  the  intercepted  snow, 
that  water  equivalent  is  merely  reduced  to 
satisfy  the  evaporation  requirement,  V^.. 
However,  if  equation  [6]  indicates  a  greater 
value  than  the  intercepted  water  equivalent,  V^, 
is  reduced  to  the  point  where  the  requirement  is 
satisfied  by  the  water  equivalent  of  the  snow, 
which  is  completely  vaporized  from  the  canopy. 

Subroutine  SNOWVAP  (fig.  4) 

During  conditions  when  the  canopy  is 
free  of  snow,  that  is,  when  time  since  the 
beginning  of  the  last  snowfall  event  is  greater 
than  2  days,  it  is  assumed  that  evaporation 
from  the  snowpack  beneath  the  trees,  v  ,  takes 
place  according  to  the  relation 


V    =  (1  -  C.)  E 
s  d  s 

when 

Cd  =  0(a  forest  opening), 
Vs  =  Es 


[7] 
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Figure  4.  —  Subroutine  SNOWVAP. 


Subroutine  EVTRAN  (fig.  5) 

Available  Soil  Water  Correction 

There  has  been  some  work  with  crops 
and  forest  ecosystems  to  indicate  that 
transpiration  decreases  as  available  water 
decreases  (Denmead  and  Shaw  1962,  Cowan 
1965,  Swanson  1967,  Swanson  1969).  Accord- 
ingly, equation  [3]  was  further  adjusted  to 
account  for  available  soil  water.  Denmead  and 
Shaw  (1962)  point  out  that,  in  porous  soils,  the 
decline  should  not  be  pronounced  until  most  of 
the  "available  water"  is  removed.  Because 
subalpine  soils  are  coarse  textured  (Retzer, 
1962),  it  was  assumed  that  transpiration  of 


dense  forest  cover  would  proceed  at  rates  given 
by  equation  [3]  until  the  soil  water  is  depleted  to 
50  percent  of  field  capacity.  Thereafter, 
transpiration  is  decreased  in  proportion  to  the 
amount  of  available  soil  water  below  one-half 
of  field  capacity.  In  open  cutover  areas,  it  was 
reasoned  that  the  absence  of  dense  vegetation 
would  enable  transpiration  to  proceed  at  rates 
given  by  equation  [3]  only  when  soil  is  at  field 
capacity.  In  the  model,  it  was  assumed  that  the 
available  soil  water  (mantle  storage)  is  5.3 
inches  in  both  the  forest  and  open,  based  on  an 
assumed  average  rooting  depth  of  4  feet,  and  a 
"wilting  point"  and  "field  capacity"  of  4  per- 
cent and  15  percent  by  volume,  respectively. 
Thus,  equation  [3]  was  expanded  to  obtain 
"actual"  evapotranspiration  in  the  forest,  E^, 
and  in  the  open,  E^^,  during  the  growing 
season  as  follows: 

For  Forest: 

E  ^  =  (1  -  R.)   (0.377M)   (E  )  [8] 
at  £  s 

For  Open: 

E      =  (1  -  R^)  (0.755M  -  3)  (E  )  [8a] 
ao  r  s 

where 

M  =  "available"  mantle  storage.  When  M  ex- 
ceeds 5.3  and  2.65  inches  in  the  open  and 
forest,  respectively,  evapotranspiration  is 
computed  by  equation  [3],  and 

Rf  =  reflectivity  of  the  forest  stand  or  open  area 
as  discussed  below. 
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Figure  5.  —  Subroutine  EVTRAN. 


Radiation  Balance  and 
Evapotranspiration  According  to 
Forest  Cover  Density 


Baumgartner  (1967)  and  Tajchman 
(1971)  have  reported  that  evapotranspiration 
from  coniferous  forests  is  greater  than  from 
open  land,  although  Tajchman  reported 
smaller  differences  between  forest  and  open 
land  than  did  Baumgartner.  Both 
Baumgartner  and  Tajchman  discussed  the 
differences  in  evapotranspiration  from  various 
cover  types  in  terms  of  the  differing  energy 
balances.  In  presenting  an  analysis  of  the 
radiation  balance  and  associated  vapor  loss, 
Baumgartner  (1967)  pointed  out  that  "the  only 
pertinent  variations  with  regard  to  the  latent 
heat  flux  are  those  associated  with 
reflectivity  ..."  Accordingly,  a  relationship 
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was  derived  between  reflectivity  and  forest 
cover  density  (Leaf  and  Brink  1973)  to  index  the 
reduction  of  evapotranspiration  as  forest  cover 
is  removed.  For  lack  of  any  field  data,  the 
following  tentative  relationships  were  assum- 
ed: 


=  0.5  - 


0.75  C 


dmx 


[9] 


where 
Cd 


the  reflectivity  of  the  forest  stand, 
natural  forest  cover  density,  expressed 
as  a  decimal,  and 

reduced  forest  cover  density,  as  a 
decimal. 


Equation  [9]  only  applies  when  <_ 


dmx 


When  > 


dmx 


,  the  reflectivity  is  given  by 


dmx 


(C,  -0.33C,  ) 
d  dmx 


[10] 


The  relationship  given  by  equations  [9]  and  [10] 
is  plotted  in  figure  6.  Note  that  when  cover 
density  Cd  =  Cc[mx»  Rf  =0.1,  whereas  when 
=0,  Rf  =  0.5.  These  values  qualitatively 
agree  with  values  given  by  Baumgartner 
(1967),  who  summarized  variations  of  absorp- 
tion coefficient  for  several  cover  types,  and 
Burroughs  (1971),  who  developed  a  shortwave 
reflectivity  model  for  lodgepole  pine  forest 
which  accounts  for  varying  stand 
characteristics  and  season  of  the  year. 


=  0.25 
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Figure  6.  —  Assumed  variation  of  reflectivity  (Rf ) 
as  a  function  of  forest  canopy  density(C^). 


Seasonal  Course  of  Transpiration 

Swanson  (1967)  observed  that  transpira- 
tion can  occur  early  in  the  snowmelt  runoff 
season  when  there  is  still  considerable  snow 
cover.  At  the  Fraser  Experimental  Forest,  dur- 
ing the  first  week  in  May  1965,  he  observed  a 
sharp  upturn  in  sap  flow  when  the  snow  cover 
still  held  an  average  of  more  than  5  inches  of 
residual  water  equivalent.  Accordingly,  a 
threshold  water  equivalent  was  assumed  in  the 
model  after  which  evapotranspiration  is  allow- 
ed to  occur.  This  threshold  is  tentatively  es- 
timated to  be  5  inches.  Evapotranspiration  is 
computed  by  equations  [8]  or  [8a]  above,  and 
varies  according  to  available  energy.  Eg, 
available  mantle  storage,  M,  and  reflectivity, 
Rf ,  all  of  which  have  been  discussed  previously. 

Once  the  evapotranspiration  re- 
quirements have  been  established  by  the  above 
adjustments,  they  are  satisfied  first  from  the 
input  and  then  from  the  soil  mantle  storage.  If 
the  requirements  would  deplete  storage  below 
the  wilting  point,  however,  all  values  are  ad- 
j\isted  to  cause  evapotranspiration  to  cease  at 
that  point. 

During  the  winter,  evapotranspiration  is 
computed  by  equations  [8]  or  [8a],  provided  the 
forest  canopy  is  free  of  snow  and  the  snowpack 
water  equivalent  is  less  than  the  critical  5 
inches.  When  the  snowpack  exceeds  5  inches, 
only  evaporation  from  intercepted  snow  and 
from  the  snow  surface  takes  place. 

Subroutine  RADBAL  (fig.  7) 

This  routine  is  essentially  identical  to 
the  routine  by  the  same  name  in  the  snowmelt 
model  (Leaf  and  Brink  1973).  The  only 
changes  necessary  were  to  modify  the 
calculations  for  year-round  processing. 
Hence,  the  modification  in  RADBAL  consists 
of  a  phase  switch,  which  indicates  optional 
methods  of  computing  the  radiation  balance. 

During  the  fall  and  winter  before  the 
diffusion  model  achieves  mathematical 
stability  (Leaf  and  Brink  1973),  only 
shortwave  radiation  is  used  to  compute 
snowmelt.  Therefore,  the  only  cold  content  in 
the  snowpack  results  from  newly  fallen  snow. 
(To  insure  snowpack  accumulation  during 
this  phase,  snowmelt  is  not  "allowed"  to  take 
place  on  days  when  the  mean  air  temperature 
is  below  0°  C).  Once  the  snowpack  depth  is 
sufficient  for  diffusion  model  stability  (4.7 
inches  of  water  equivalent),  the  phase  switch 
is  reset  to  compute  snowmelt  according  to  the 
snowmelt  simulation  model  (Leaf  and  Brink 
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1973).  In  other  words,  when  the  snowpack 
reaches  the  specified  depth  in  the  early 
winter,  the  switch  is  set  to  "melt  when  ready." 
The  snowpack  will  continue  to  accumulate  in 
this  phase  until  the  normal  all- wave  radiation 
balance  produces  snowmelt.  The  switch  is 
reset  to  the  snow  accumulation  phase  at  the 
end  of  the  snowmelt  season  (when  snowpack 
water  equivalent  is  zero). 


RADBflL  J 


COMPUTE  SW 

RADIATION 

COMPONENT 


ACCOM. 


(SUMMER, FALL 


BETWEEN 


WINTER) 


SET  PHASE  TO 
"IN  BETWEEN" 
AND  INITIALIZE 
DIFFUSION  MODEL 


C 


MELT 
(WINTER, 
SPRING) 


YES 


RETURN 


> 


SET  PHASE 
INDICATOR  TO 
"MELT  WHEN 
READY" 


CHECK  WITH 
LINKING 
ROUTINE 


Figure  7.  —  Subroutine  RADBAL. 


Subroutine  SNOWED  (fig.  8) 

The  only  difference  between  this  routine 
and  its  counterpart  in  the  snowmelt  model 
(Leaf  and  Brink  1973)  is  the  inclusion  of 
interception  calculations.  The  amount  of 
snow  intercepted  on  a  given  day  is  computed 
as  a  percentage,  which  is  determined  by  the 
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Figure  8.  —  Subroutine  SNOWED. 


vegetation  type  as  discussed  in  subroutine 
CANVAP.  The  total  amount  which  may  re- 
main on  the  canopy  after  several  consecutive 
snow  events  is  also  determined  by  vegetation 
type. 


Peripheral  Routines 

The  peripheral  routines  are  not  flow- 
charted  or  discussed  in  detail  on  an  individual 
basis,  since  they  are  primarily  utility  routines 
for  input,  output,  and  maintenance  of  con- 
tinuity between  simulation  intervals.  Sub- 
routines GENDATA,  RADCOMP,  and 
RDMSTR  are  all  concerned  with  input,  and 
use  various  analyses  to  generate  the  daily 
input  for  each  substation  from  observed  base 
station  data.  Subroutines  ETCODE  and 
PLOTTER  are  output  routines  which  aid  in 
the  interpretation  of  results. 

All  of  the  above  routines  are  utilized  on 
each  run  of  the  model,  but  to  save  time  and 
core,  most  output  options  are  included  as 
overlays,  only  one  of  which  is  selected  at  a 
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time  to  occupy  core  and  process  the  data. 
Each  overlay  consists  of  a  control  routine,  a 
normal  simulation  routine,  a  write  routine, 
and  any  watershed  management  alternative 
simulation  routines  that  may  be  needed. 

A  complete  listing  of  the  model  describ- 
ed in  this  report  is  included  in  the  appendix. 


Applications 

We  have  used  the  hydrologic  model  to 
simulate  area  snowmelt  and  water  yield  from 
the  667-acre  Deadhorse  Creek  watershed  at 
the  Fraser  Experimental  Forest  (Leaf  1971; 
Leaf  and  Brink  1972,  1973).  Physical 
characteristics  of  this  watershed  vary  from 
low-elevation  (9,300  ft.  m.s.l.)  south  slopes  in 
lodgepole  pine  forest  to  high-elevation  (11,000 
ft.)  north  slopes  in  spruce-fir.  Simulations 
from  10  subunits  on  the  basin  were  weighted 
according  to  the  percentage  of  the  total  area 
each  represents  to  generate  the  watershed 
composite.  Each  subunit  was  selected  accord- 
ing to  forest  cover  type  and  density,  slope, 


aspect,  and  average  elevation.  Simulations  of 
the  1963-70  water  years  (October  1  to 
September  30)  indicate  to  us  that  the  model 
represents  the  inherent  hydrologic 
characteristics  of  Colorado  subalpine 
watersheds.  Sample  output  in  the  form  of  10- 
day  summations  for  the  1967  water  year  is 
shown  in  table  1.  Figure  9  summarizes  10-day 
fluctuations  of  several  hydrologic  variables 
for  the  1965  water  year.  Figure  10  compares 
simulated  and  observed  annual  water  yields 
from  Deadhorse  Creek  for  the  1964-71  record 
period. 

We  have  predicted  the  change  in  rate 
and  seasonal  time  distribution  of  snowmelt 
resulting  from  clearcutting  small  openings  in 
old-growth  forest  with  the  snowmelt  portion 
of  the  model  (Leaf  and  Brink  1972).  With  the 
increased  water  balance  simulation  capabili- 
ty, we  plan  to  develop  the  model  into  a  useful 
tool  for  predicting  the  hydrologic  conse- 
quences of  several  resource  management 
practices.  In  the  Colorado  subalpine  zone, 
these  include  weather  modification  and 
timber  harvesting. 


Figure  9.  —  Simulated  10-day  fluctuations  of  several  hydrologic  components  during 
the  1965  water  year  on  Deadhorse  Creek,  Fraser  Experimental  Forest. 
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Figure  10.  —  Simulated  versus  observed  annual 
runoff  on  Deadhorse  Creek,  1964-71. 
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Appendix  I:  Routines  for  Water  Balance  Model 


Program  WBMODEL 


DVEBL«y  (OLAYS.0,0) 

PROGRAM  HBMODEL   I  I NPUT , OUTPUT, PLOTS .TAPESi I NPUT ,T»PE6=0UTPUT , 
1  TAPEll-PLOTS) 

C  THIS   IS  THE  CONTROLLING   ROUTINE   FOR   THE   HATER  BALANCE   MODEL.  THE 

C   MODEL   IS  OVERLAYED  TO  SAVE  TIME  AND  MEMORY.  THE  CORE  OF  THE  MODEL 

C   IS  SUBROUTINE  HATBAL  AND   ITS  RELATED  ROUTINES.   ALL  OF  THESE,  PLUS 

C   SEVERAL   I/O  ROUTINES  WHICH  HAVE  NO  OPTIONS,   ARE   INCLUDED   IN  THIS 

C   MAIN  OVERLAY   SO  THEY  ARE   AVAILABLE   TO  THE   OVERLAYS  WHICH  ARE 

C   SELECTED  BY  OPTION.     THE  FIRST  OVERLAY  LOADED  ESTABLISHES  THE 

C   WATERSHED  DESCRIPTORS  AND  PARAMETERS  AND  INDICATES  THE  OVERLAYS 

t  "TO  BE   SELECTED  ACCORDING  TO  THE   OUTPUT   OPTIONS.      THE   NEXT  OVERLAY 

C   TO  BE  LOADED  WILL  CONTAIN  THE  OUTPUT  ROUTINE  SELECTED  AND  ITS 

C   NORMAL    SIMULATION  ROUTINE.      IT   WILL   THEN   LOAD  ONE   OF  ITS 

t   SECONDARY  OVERLAYS  WHICH  HILL  CONTAIN   THE   MAIN  OPERATING  PROGRAM 

C   FOR   THIS   RUN  AND   ANY   ALTERNATIVES   THAT   ARE   TO   BE  INCLUDED. 

C  

C  THE  CORE  OF  THE  MODEL   (SUBROUTINE  WATBAL,   ET  AL I    IS  DESIGNED  TO 

C   FUNCTION  AS  AN   INDEPENDENT  UNIT,   TOTALLY   UNCONCERNED  WITH   I/O  AND 

C   MAINTENANCE^   OF  CONTINUOUS  OR  STATIC  CONDITIONS.     THE  ONLY  LINKS 

C   WITH  THE  OPERATING   PROGRAM,    ITS  ALTERNATIVES  AND   I/O  ROUTINES, 

C   ARE   THE  FORMAL   PARAMETERS   AND  COMMON  BLOCK   /WATRBAL/.  EACH 

C   ROUTINE   WHICH  UTILIZES   THE   CORE  MUST   MAINTAIN   ITS   OWN   SET  OF 

C   CONTINUOUS  CONDITIONS   AND  MAKE   THEM   AVAILABLE   AT   THE   PROPER  TIME 

C   FOR  USE  BY  THE  CORE.     THE  BLANK  COMMON  AND  LABELLED  COMMON  BLOCKS 

C   ESTABLISHED   HERE    IN    THE    MAIN   OVERLAY    ARE    PRIMARILY    FOR    THE    USE  OF 

(;  The   NORMAL   SIMULATION  ROUTINE.     THEREFORE,    SIMILAR   LOCATIONS  MUST 

C   BE  SET  ASIDE  FOR   MAINTENANCE     BY  ANY  AND  ALL  ALTERNATIVES. 


COMMON  AIRTEmCI25,6I 
COMMON  CALDEF (25) ,COMAXI25l ,CDVDEN(25 I 
COMMON  DREADy(25) 

COMMON  ENGBAL(25I,ET,ET0AILY|25,12) 

COMMON  FREEWAT(25) 

COMMON   I SOTHRM(25,20l 

COMMON  LASTUSD(251 ,LEVELl,LEvEL2 

COMMON  MMDD 

COMMON  NDAYSN0( 25) , NDI VSBL , NSUB ,NYE ARS 
COMMON  0NTREES(25) 

COMMON  PEAKPPT(25,20),PEAKWE(  25,20) , PHASE! 25 ) , POTENT( 24  I , 
1       .  PREWE0V(25) 
COMMON  RECHRG(25) 

COMMON  SIM  TEMK  25. 3), SUB  10(25, 51  ,  SLPASPI  25,  2<,) 

COMMON  TC0EFFI25) ,THRSHL0( 25) ,TOPLOT( 11 ) 

COMMON  VEGTYPE(25) 

COMMON   WEIGHTI25I ,WSHE0IDI6) 

COMMON  YEARS(20I ,YYHMDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TCTPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDID 

INTEGER   YEARS, YYMMDD 

CDMHON/CMPSIT0/COHPSI16),YRTOT(5) 

COMMON/MASTER/ DATE ( 3 ) , TMXM STR , TMNMSTR , PPTMSTR , PPTONOW , OB SHYDR , 
1  POTRAD,HSTREOF, I YR 
INTEGER  DATE 

COMH0N/WATRBAL/ETFR0M,EVAP0TR,GENR0,PRECIP,RADIN,RADLWN,RADSWN, 
1  TEHPMAX,TEMPHIN,WATERIN 
DATA  C0MPS,YRT0T/21«0.0/ 
DATA  IYR,HSTRE0F/1 ,0/ 
CALL  OVERLAY  (5h0LAYS,1 ,0) 
CALL  OVERLAY  ( 5H0LAYS, LEVEL  1 , 0 ) 

C  FILE  -PLOTS-   IS  COPIED   To  -OUTPUT-   BY   A  STANDARD  MONITOR  COPY 

END 

Subroutine  WATBAL 


SUBROUTINE  WATBAL  (FI,F2,F3,I1,F4,F5,I2,I3,F6,U.F7,F8,F9,F10,FU, 
1  F12,F13,I5I 

C  THIS  SUBROUTINE   IS  THE  MAIN  ROUTINE  OF  THE  WATER  BALANCE  MODEL.  IT 

C   RECEIVES  THE   DRIVING,    STATIC   AND  CONTINUOUS  VARIABLES   FROM  THE 

C   OPERATING  ROUTINES,  CONTROLS  THE  COMPUTATIONS  ON  THEM,  AND 

C   RETURNS   THE   NEW  VALUES   FOR  THE  CONTINUOUS   VARIABLES  AND  THE 

C   RESULTS  OF  THIS   INTERVAL.     SEE  THE  REPLACEMENT  STATEMENTS  BELOW 

C   FOR  THE  VARIABLE  DEFINITIONS  OF  THE  PARAMETERS 

C 

C  DICTIONARY  OF  WATER  BALANCE  COMMON  BLOCKS 

C 

C  AVETEHC  -  THE  MEAN  TEMPERATURE  FOR  THE   INTERVAL   IN  DEGREES  C 

C  BASTEMF  -  BASE  TEMPERATURE  DEGREES  FARENHEIT,   RAIN  TURNS  TO  SNOW 

C  CALDEF  -  THE  CALORIE  DEFICIT   IS  THE  NUMBER  OF  CALORIES  NEEDED 

C  TO  BRING  THE   SNOWPACK   TEMPERATURE   TO   ZERO  DEGREES 

C  CENTIGRADE   (NOTE   SHOULD  BE  MADE  THAT   IT   IS  A  POSITIVE 

C  QUANTITY) 

C  COVDEN  -  THE  COVER  DENSITY  IS  THE  FRACTION  OF  THE  GROUND  OR  SNOW 

C  SURFACE  SHADED  FROM  DIRECT  SUNLIGHT  OR  RADIATION 

C  DREADY  =  0,  DIFFUSION  MODEL   (SUBROUTINE  DIFMOD)   NOT  INITIALIZED 

C  =1,   DIFFUSION  MODEL   INITIALIZED  AND  READY   FOR  SNOWPACK 

C  TEMPERATURE  SIMULATION 

C  =  -I,  DIFFUSION  MODEL  HAY  NOT  BE  USED 

C  ENGBAL   -  THE   TOTAL  CALORIC   INPUT   TO  OR  LOSS   FROM  THE  SNOWPACK 

C  DURING   AN   INTERVAL.      IT    IS  THE   ALGEBRAIC   SUM   OF  THE 

C  ENERGY   INVOLVED  WITH  THE   PRECIPITATION  AND   THAT  OF 

C  THE   RADIATION  BALANCE 

C  ENGBALl  -  THE   VALUE   OF   -ENGBAL-  AT   THE   END  OF  THE  LAST  INTERVAL 

C  ETFROM  =  I,   EVAPORATION   IS   FROM  THE  CANOPY 

C  =2,  EVAPORATION  IS  FROM  THE  SURFACE  OF  THE  SNOWPACK 

C  '         EVAPOTRANSPIRATION   is  from   SNOWMELT,   rain  or  THE 

c  soil  mantle  storage 


C  EVaPOTR   -   WHEN  FIRST   RECEIVED,   THIS   VARIABLE   IS  THE  POTENTIAL 

C  EVAPOTRANSPIRATION    AS   COMPUTED    BY   THE    HAMON  METHOD 

C  AND   ADJUSTED  FOR   AVAILABLE   RADIATION.      AFTER  ACTION 

C  IS   TAKEN  BY  THE   WATER   BALANCE   ROUTINES,   THE  ORIGIMAL 

C  VALUE   HAS   BEEN  ADJUSTED  FURTHER   BY   THE  METHODS 

C  DISCUSSED   IN   SUBROUTINES  CANVAP,   EVTRAN,   AND  SNOWVAP. 

C  IT   THEN  REPRESENTS   THE   EVAPOTRANSPIRATION  DURING  THIS 

C  INTERVAL 

C  FREEWAT  -   THE   FREE   WATER  BEING  HELD  BY   THE  SNOWPACK 

C  LASTUSD  -   AN   INDICATOR   USED   IN   FUNCTION  PACKREF   TO  DETERMINE 

C  WHICH   REFLECTIVITY  FUNCTION  TO  USE 

C  NDAYSNO   -    THE    NUMBER   OF    DAYS    SINCE    NEW    SNOW    HAS  FALLEN 

C  ONTREES  -   THE   VOLUME   OF    INTERCEPTED   SNOW  REMAINING  ON   THE  CANOPY 

C  PHASE   =   -1.    THE   SNOWPACK  HAS   BARELY  ACCUMULATED  TO  THE  POINT  WHERE 

C  THE   DIFFUSION  MODEL   IS   STABLE   AND  MAY   BE  USED  TO 

C  CONTROL    THE    SNOWPACK  TEMPERATURE 

C  =0,   THE  SNOWPACK   IS  ACCUMULATING  AND  HAS  NOT  YET  REACHED 

C  A  DEPTH  WHICH  WILL   PROVIDE  STABILITY  FOR  THE 

C  DIFFUSION  MODEL 

C  -   I,    THE    SNOWPACK   HAS   REACHED  A   SUFFICIENT   DEPTH  TO  ALLOW 

C  THE   DIFFUSION  MODEL   TO  CONTROL   THE   PACK  TEMPERATURE 

C  UNTIL   THE   MELT   SEASON,   WHEN  THE   RADIATION  ROUTINES 

C  RESUME  CONTROL  TO  GOVERN  THE  MELT  PHASE 

C  PRECIP  -  OBSERVED  PRECIPITATION   IN  INCHES 

C  PREWEQV   -  PREDICTED   WATER  EQUIVALENT   OF   THE   SNOWPACK   IN  INCHES 

C  RADIN  -   RADIATION   IN   IS   THE   TOTAL   INCIDENT   SHORT  WAVE  RADIATION 

C  RADLWN  -   NET   LONG   WAVE   RADIATION   IS   THE   ALGEBRAIC   SUM  OF   THE  LONG 

C  WAVE   RADIATION   FROM  THE   FOREST   AND  THE   LONG  WAVE 

C  RADIATION  LOST  BY   THE   SNOWPACK   TO  THE  CANOPY 

C  RADSWN  -   THE   CALORIC   INPUT   TO  THE   PACK   BY   THE  NET   SHORT  WAVE 

C  RADIATION 

C  RECHRG  -  THE  RECHARGE  REQUIREMENTS,  OR  SOIL  MANTLE  STORAGE  DEFICIT 

C  SIMTEMl  -  AN   ARRAY  USED  PRIMARILY   IN   SUBROUTINE   DIFMOD   IN  THE 

C  SIMULATION  OF  THE  AVERAGE  SNOWPACK  TEMPERATURE. 

C  TO   INSURE   STABILITY   OF   THE  DIFFUSION  MODEL,  THE 

C  DAY    IS   PARTITIONED   INTO   12  HOUR   INTERVALS,  AS 

C  DISCUSSED   IN   SUBROUTINE   DIFMOD.     THIS   ARRAY  STORES 

C  THE   CONDITIONS   PRESENT   DURING  THIS    INTERVAL   FOR  USE 

C  IN  THE   SIMULATION  ON   THE  NEXT    INTERVAL.      LOCATION  1 

C  STORES  THE  AVERAGE   AIR  TEMPERATURE   (ASSUMED  TO  BE 

C  The   SURFACE  TEMPERATURE  OF  THE  SNOWPACK),  LOCATION  2 

C  IS   THE   SNOWPACK   TEMPERATURE   AT   A  NODE  MIDWAY 

C  BETWEEN  THE  SURFACE  AND  THE  GROUND,   AND  LOCATION  3  IS 

C  THE   GROUND  TEMPERATURE. 

C  TCOEFF  -  THE   TRAN SMI  SI  V I T Y  COEFFICIENT  USED  TO  ESTIMATE  THE  NET 

C  SHORT  WAVE   RADIATION  REACHING  THE  SNOWPACK.  SEE 

C  REIFSNYDER   AND  LULL,   RADIANT   ENERGY   IN  RELATION  TO 

C  FORESTS,   USFS  TFCH.   BUL   13*<i.  1965. 

C  TEMPMAX  -   THE   MAXIMUM   TEMPERATURE   DURING  THE    INTERVAL    IN  DEGREES 

C  FARENHEIT  •  ' 

C  TEMPMIN  -  THE  MINIMUM  TEMPERATURE  DURING  THE   INTERVAL   IN  DEGREES 

C  FARENHEIT 

C  THRSHLD  -  THE   THRESHOLD  TEMPERATURE  FOR  DETERMINING  WHETHER  OR  NOT 

C  TO  RE-INITIALIZE   THE   REFLECTIVITY  FUNCTION  WHEN 

C  THERE   IS  A   SNOW  EVENT.     IF  THE  MAXIMUM  TEMPERATURE  IS 

C  GREATER  THAN  THE  THRESHOLD  VALUE  DO  NOT  RE- IN  IT  I AL IZE 

C  THE   FUNCTION  REGARDLESS   OF  THE  PRECIPITATION 

C  WATERIN  -   THE   SUM  OF   ANY   SNOWMELT   AND  ANY  RAIN   WHICH  PROVIDES 

C  DIRECT   INPUT  TO  THE  WATER  BALANCE 

C 

COHMON/ONLYCOR/   A VE TE MC , BASTEMF , CALDEF, COM AX, COVDEN, DREADY, ENGBAL , 


1  ENGBALl, FREEWAT, LASTUSD, NDAYSNO, ONTREES, PHASE, PREWEQV, RECHRG, 

2  SIHTEMU  3)  ,  SI  MTEM3,  TCOEFF  ,  THRSHLD,  VEGTYPE 
INTEGER  DREADY, PHASF, VEGTYPE 

C OMMON/ W AT R BAL /E TF ROM, EVAPOTR,GENRO, PRECIP, RAD  IN, RADLWN. RADSWN, 
1    TEMPMAX. TEMPMIN. WATERIN 
DATA   A vETEMC. BASTEMF .CALDEF .CDMAX, COVDEN, DREADY, ENGBAL, ENGBALl, 

1  FREE  WAT. LASTUSD,NOAYSNO,ONTREES, PHASE, PREWEQV, RECHRG, SIMTEMl, 

2  SI MTEM3,TCOEFF,THRSHLD,VEGTYPE/0. 0,35.0, 3»0.0, 0, 2»- 1 .0,0.0, 2*0, 

3  0.0,0.6*0.0.1.0.0.0.1/ 


C  OBTAIN   THE   STATION  DESCRIPTORS 

COVDEN  =  F3 
CDMAX  =  F2 
TCOEFF  =  F12 
VEGTYPE  =  15 

C  RECALL  THE  CONTINUOUS  VARIABLES  NECESSARY  FOR  THE  OPERATION  OF  THE 

C   MODEL   DURING   THIS  INTERVAL 


CALDEF   =  Fl 
DREADY  =  II 
ENGBALl  =  F* 
FREEWAT  =  F5 
LASTUSD  =  12 
NDAYSNO   =   13  t  1 
ONTREES   =  F6 
PHASE   =  I* 
PREWEQV   =  F7 
RECHRG  =  F8 
THRSHLD  =  F13 


IF(DREADY)  20.20.10 
10   SIMTEMK  1)    =  F9 
SIMTEMl (2)    =  FIO 
SIMTEM1(3)    =  Fll 

C  AVETEMC  =   ( ( ( TEMPNAX-32 )♦( TEHPMIN-32) l/2l»l5/9) 

20  AVETEHC   =   (TEMPMAX   t   TEMPMIN  -  b'l.O)    •  0.2777777778 

C  START   THE   ENERGY  BALANCE  AND  THE   INPUT   AT   ZERO  FOR   THIS  INTERVAL 

ENGBAL   =  0.0 
WATERIN   =  0.0 

C  IF   THERE   IS  NO  PRECIP.   THERE   IS  NO  NEED  TO  PASS  THROUGH  THE 

C   CLASSIFICATION  STATEMENTS 

IF(PRECIP)  90.90.30 
C  SEE   IF    THE   PRECIP   IS  ALL  SNOW 
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IF  I  TEHPM1N.L6. 32.0. OR. T6MPH4x.lt. 6ASTEHFI   GO  TO  60 

C  SEE   IF  ANY  OF   IT   IS  SNOW 

IF(TEHPH|N  -  BAST6MF)  10,50,50 
<iO  CALL  MIXTURE 
GO  TO  90 

C  THIS    IS   A  RAIN  EVENT.      IF   THERE    IS  NO  PACK.    THE  RAIN   IS  DIRECT 

C   INPUT  TO  THE  HATER  BALANCE.     BUT   IF  THERE   IS  A  PACK,  DETERMINE 

C   THE   EFFECTS  OF   THE  RAIN 

50   IF(PREHEQV]  60,60,70 

60  HATERIN  =  PRECIP 
GO  TO  160 

70  CALL  Rained  iavetemcprecip) 

GO   TO  90 
C  THIS    IS    A   SNOW  EVENT 

BO  CALL  Snowed  (amini  (avetehc.o.O) ,precip) 

C  IF   THERE   IS  SNOW  ON  THE  TREES,   EVAPORATE  ONLY  FROM  THE  CANOPY 

90  IF(ONTREES)  130,130,100 
100  CALL  CANVAP 

C  ON  THE  FIRST  DAY  AFTER  FRESH  SNOW,   ASSUME  TURBULENCE  HAS  REMOVED 

C   ANY  REMAINING  INTERCEPTED  SNOW  AND  ADDED  IT  TO  THE  PACK 

IFINDAYSNO  -   11  120,110,110 
110  PREWEOv  =  PREhEOV  *  ONTREES 
ONTREES  =  0.0 

C  IF   THERE    IS   NO  SNOWPACK,   BYPASS   THE  RADIATION  ROUTINES 

120   IF(PREWEOV)  190,190,180 

C  DETERMINE  WHETHER  TO  SATISFY  THE  EVAPOTRANSPIRATION  REOUIREMENTS 

C   UNDER  GROWING  SEASON  OR  WINTER  CONDITIONS 

130   IF(PREHEQV)  160,160,140 

140   IFIPREWEQV  -  5.0)  150,150,170 
C  USE  THE  GROWING  SEASON  ROUTINES  TO  INCLUDE  TRANSPIRATION 

150  CALL  RADBAL 

C  ADD  -WATERIN-  TO  THE  RECHARGE  REQUIREMENTS  SO  THE  ET  ROUTINE  CAN 

C   OPERATE  ON  THE   INPUT  AS  WELL  AS  THE  STORAGE 

160  RECHRG  =  RECHRG  ♦  WATERIN 

CALL  EVTRAN 

GO  TO  200 

C  USE  THE  WINTER  ROUTINES  TO  EVAPORATE  FROM  THE  SNOWPACK  SURFACE 

170  CALL  SNOWVAP 

180  CALL  RADBAL 
C  ADD  -WATERIN-  TO  THE  RECHARGE  REQUIREMENTS 

190  RECHRG  =  RECHRG  ♦  WATERIN 

C  IF   THE  RECHARGE  REOUIREMENTS  WERE  SATISFIED,   THE  EXCESS  IS 

C   CONSIDERED  TO  BE  GENERATED  RUNOFF 

200   IFIRECHRGI  220,220,210 

210  GENRO  =  RECHRG 
F8  =  0.0 
GO  TO  230 


220 

GENRO  =  0.0 

F8  = 

RECHRG 

230 

II  = 

OREADY 

Fl  = 

CALDEF 

F4  = 

ENGBAL 

F5  =■ 

FREEHAT 

12  = 

LASTUSO 

13  » 

NDAYSNO 

F6  > 

ONTREES 

F7  = 

PREWEQV 

C  WHEN  THE  PACK  IS  GONE,   RESET  THE  PHASE  INDICATOR 

IF(PREWEQV)  2'>0,2'.0,250 
240  14  =  0 

RETURN 
250  14  =  PHASE 

IFIDREADYI  270,270,260 
260  F9  =   SIMTEMll I ) 

FIO  =  SIMTEM1I21 

Fll   =  S1MTEM1(3I 
270  RETURN 

END 


Subroutine  CALIN 

SUBROUTINE  CALIN  (CALORINI 

C  THIS  SUBROUTINE  COMPUTES  THE  EFFECTS  OF  THE  CALORIC   INPUT  ON  THE 

C   SNOWPACK 

C  HOLOCAP  -  THE  FREE  WATER  HOLDING  CAPACITY  OF  THE  SNOWPACK 

C  (ASSUMED  TO  BE  FOUR  PERCENT  OF  THE  WATER  EQUIVALENT! 

C  

COMMON /ONLYCOR/  AVE TEMC , BAST EMF , CALDEF , COM  AX , COVOEN , DR EADY , ENGBAL 

1  ENGBAL I, FREE WAT, LASTUSO, NDAYSNO, ONTREES, PHASE, PREHEOV, RECHRG, 

2  SIMTEMK  31 ,SIMTEM3,TC0EFF,THRSHL0, VEGTYPE 
INTEGER  DREADY, PHASE .VEGTYPE 

C OHMON /WATRBAL/ETFROM,EVAPOTR, GENRO, PRECIP, RAD  I N,RADLWN,RAOSHN, 
1   TEMPMAX,TEMPHIN, WATERIN 
C  ADD  THESE  CALORIES  INTO  THE  ENERGY  BALANCE 

ENGBAL   =  ENGBAL   ♦  CALORIN 
C  SEE    IF   A  CALORIE   DEFICIT   EXISTS    IN  THE  PACK 

COMPARE   =  CALORIN  -  CALDEF 

IF(COMPARE)  10,20,30 

C  THERE   IS  A  CALORIE  DEFICIT,   BUT  THE   INPUT  DID  NOT  COMPLETELY 

C   WIPE    IT  OUT.      ALL   OTHER  CONDITIONS   ARE  UNCHANGED 

10  CALDEF  =  -  COMPARE 

RETURN 

C  THE  CALORIE  DEFICIT  HAS  WIPED  OUT,   BUT  ALL  OTHER  CONDITIONS  ARE 

C   UNCHANGED 

20  CALDEF  =  0.0 
RETURN 

C  ANY  DEFICIT  WHICH  DID  EXIST  WAS   WIPED  OUT.      COMPUTE   THE  POTENTIAL 

C   MELT  FROM  THE  REMAINING  CALORIES   ( CALOR I ES/ ( 80 .0  *  2.54)1 

30  POTMELT  =  COMPARE/203.2 
CALDEF   =  0.0 

C  IF   THE    INPUT  WAS   ENOUGH  TO  MELT   THE   WHOLE  PACK,   CONTRIBUTE  THE 

C   WATER   EQUIVALENT   TO  THE   SNOWMELT   AND   ZERO   ALL  CONDITIONS 

IFIPOTHELT.lt. PREHEOV-FREEHAT)   GO  TO  40 

WATERIN   =   WATERIN  ♦  PREWEQV 

PREHEQV  =  0.0 


FREEHAT   =  0.0 
RETURN 

C  DEPLETE   THE   ICE  PACK  BY  THE  AMOUNT  MELTED  AND  CONTRIBUTE  THAT 

C   AMOUNT  TO  THE   FREE  HATER 

40  FREEWAT  =  FREEWAT  ♦  POTMELT 

C  COMPUTE   THE   NEW  HOLDING  CAPACITY  OF   THE   PACK  AND  COMPARE   IT  WITH 

C   THE  FREE  WATER  TO  SEE   IF  SNOWMELT   IS  PRODUCED 

HOLDCAP  =  0.04  •   (PREWEQV  -  FREEWAT) 

COMPARE  =  FREEWAT  -  HOLDCAP 

IF(COMPARE.LE.O.O)  RETURN 

C  THE  SNOWMELT  CONTRIBUTED  IS   IN  -COMPARE-.     REDUCE  THE  FREE  HATER 

C   TO  LEAVE  A  PRIMED  PACK  AND  REDUCE  THE  PREDICTED  WATER  EQUIVALENT 

PREHEQV  =  PREWEQV  -  COMPARE 

WATERIN  =  WATERIN  t  COMPARE 

FREEHAT  =  HOLOCAP 

RETURN 

END 


Subroutine  CALOSS 

SUBROUTINE   CALOSS  (CALOUTt 

C  THIS  SUBROUTINE  COMPUTES  THE  EFFECTS  OF  THE  CALORIC  LOSS  ON  THE 

C   SNOHPACK 

COMMON /ONLYCOR/  AVE TEMC, BASTEMF , CALDEF, COHAX, COVOEN, OREADY , ENGBAL, 

1  ENGBAL I, FREE WAT, LASTUSO, NDAYSNO, ONTREES, PHASE, PREHEQV, RECHRG, 

2  SIHTEMl (3), SI HTEH3,TC0EFF,THRSHLD, VEGTYPE 
INTEGER   OREADY, PHASE , VEGTYPE 

COMHON/WATRBAL/ETF ROM, EVAPOTR, GENRO, PRECIP, RAD  IN, RADLWN.RADSHN, 
1   TEMPMAX.TEHPHIN, WATERIN 
C  ADD  ALGEBRAICALLY  THESE  CALORIES  INTO  THE  ENERGY  BALANCE 

ENGBAL  =  ENGBAL  +  CALOUT 

C  SEE   IF   THERE   IS  ANY  FREE  WATER  IN  THE  PACK.     IF  NOT,   THE  LOSS  IS 

C   JUST  CONTRIBUTED  TO  THE  CALORIC  DEFICIT  OF  THE  SNOHPACK. 

C   REHEMBER  THAT  -CALOUT-  IS  NEGATIVE 

IFIFREEHAT.GT.0^0)   GO  TO  10 

CALDEF  =  CALDEF  -  CALOUT 

RETURN 

C  COMPUTE  THE  CALORIC  LOSS  NECESSARY  TO  FREEZE  ALL  OF  THE  FREE  HATER 

C   (FREE   HATER   •   80.0  »  2.54) 

10  CALNEEO  =  FREEWAT  *  203.2 

C  NOW  COMPARE  THAT  NECESSARY  LOSS  WITH  THE  ACTUAL  LOSS.     IF  THEY  ARE 

C   THE  SAME,    THE  FREE  WATER   IS  WIPED  OUT  BUT  NO  OTHER  CONDITIONS  ARE 

C   ALTERED 

COMPARE  =  CALOUT  ♦  CALNEED 
IF(COHPARE)  20,30,40 

C  THE  LOSS  WAS  MORE   THAN  ENOUGH  TO  FREEZE  IT.     THE  BALANCE  CREATES 

C   AN  ENERGY  DEFICIT   IN  THE  PACK  AND  THE  FREE  HATER   IS  WIPED  OUT 

20  CALDEF  =  -  COMPARE 
30  FREEHAT  =  0.0 
RETURN 

C  ONLY  PART  OF  THE  FREE  WATER  FROZE.     COMPUTE  THE  BALANCE  REMAINING 

C   BALANCE  =  EXISTING  FREE  HATER  -  AMOUNT  FROZEN,  WHERE 

C   AMOUNT  FROZEN  =  CALOR I ES/ ( 80. 0  •  2.54) 

40  FREEHAT  =  FREEHAT  ♦   ( CALOUT/ 203 . 2 1 

RETURN 

END 


Subroutine  CANVAP 

SUBROUTINE  CANVAP 

C  COMPUTE  THE  EVAPORATION  FROM  THE  INTERCEPTED  SNOW  AS  A  FUNCTION  OF 

C   THE  CANOPY  COVER  DENSITY 

COMHON /ONLYCOR/  A VETEHC , BASTEMF , CALDEF , COM AX, COVOEN, OREADY, ENGBAL . 

1  ENGBAL 1. FREEHAT, LASTUSO, NDAYSNO, ONTREES, PHASE, PREHEQV, RECHRG, 

2  SIMTEMl(3) ,SIHTEM3,TC0EFF,THRSHLD, VEGTYPE 
INTEGER  OREADY, PHASE, VEGTYPE 

COMHON /WATRBAL/ETFROH, EVAPOTR, GENRO, PRECIP, RAD  I N,RADLWN,RADSWN, 
1   TEHPMAX,TEMPHIN, WATERIN 
ETFROH  =  1.0 

EVAPOTR  -  EVAPOTR/COVDEN 

ONTREES  =  ONTREES  -  EVAPOTR 

IFIONTREESI  10,20,20 
10  EVAPOTR   =  ONTREES  ♦  EVAPOTR 

ONTREES  =  0.0 
20  RETURN 

END 

Subroutine  DIFMOD 

SUBROUTINE  DIFMOD 

C  THIS  SUBROUTINE  HAS  DERIVED  FROM  PROGRAM  SIMTEM,   A  SNOWPACK 

C   TEMPERATURE  DIFFUSION  MODEL  DEVELOPED  BY  LEAF  (1970  STUDY  PLAN 

C   FS-RM-1602,  NO.   224,   RMF*RES).     USING  THE  AVERAGE   SURFACE  TEMP 

C   AND  THE  GROUND  TEMP  AS  BOUNDARY  CONDITIONS,  THE  NEW  AVERAGE 

C   SNOHPACK   TEMPERATURE   IS  CALCULATED 

COHHON/ONLYCOR/   AVE TEMC , BASTEMF , CALDEF , COM AX , COVOEN , ORE ADY , ENGBAL , 

1  ENGBALl, FREE WAT, LASTUSO, NDAYSNO, ONTREES, PHASE, PREWEQV, RECHRG, 

2  SIMTEHl(3 ) ,SIMTEH3,TC0EFF,THRSHL0,VEGTYPE 
INTEGER  DREADY, PHASE, VEGTYPE 

COMMON/ WAT RBAL/ETF ROM, EVAPOTR, GENRO, PRECIP, RAD  I N,RADLWN,RADSWN, 
1   TEHPMAX,TEHPMIN, WATERIN 

C  

C  DICTIONARY 

C 

C  CONSTl  -  THE  FIRST  CONSTANT   IN  THE  EQUATION  FOR  THE  SIMULATION 

C  C0NST2  -  THE  SECOND  CONSTANT  IN  THE  EQUATION  FOR  THE  SIMULATION 

C  H  -  THE  DISTANCE  BETWEEN  NODES   (CORRESPONDS  TO  THE  -H-   IN  THE 

C  STUDY  PLAN) 

C  

C  COMPUTE  THE   DENSITY  OF   THE   SNOWPACK   (THE  FUNCTION  WAS  DERIVED  FROM 

C   OBSERVED  CONDITIONS  ON  THE  FRASER  EXPERIMENTAL  FOREST) 

DENSITY  =   IEXP((0.0179  *  PREHEQV)   ♦  3.02))/l00.0 
C  COMPUTE   The  distance  between  the  two  nodes   in  CENTIMETERS 
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C  DEPTH  =  PREWEOV/DENSITy 

C  H  =  (0EPTH/2)»2.5'. 

H  =   IPREWEQV/DENSITrl   »  1.27 

C  THE   THERMAL  OIFFUSIWITY   IS  CALCULATED  FROM   THE  FUNCTION 

C   KV   =   0.01/112.751   -  DENSITY)'  0.48).     MATHEMATICAL  STABILITY 

C   REQUIRES  THAT   THE   VALUE   OF   THE   QUANTITY   (INTERVAL    IN   SECONDS  « 

C  KV/H*»2)   BE   LESS   THAN  0.5.     WHEN  A  2<t   HOUR    INTERVAL    IS   USEO,  THE 

C   SNOW  DEPTH  MUST  EXCEED  30   INCHES   120  PERCENT  DENSITY!   To  ACHIEVE 

C   STABILITY.      In  ORDER   TO   INSURE   STABILITY   WITH  SOMEWHAT  SHALLOWER 

C   PACKS    (ABOUT   18    INCHES).    THE   DAY   IS  DIVIDED   INTO  2   TIME  INTERVALS 

C  OF   12  HOURS   (<.3200  SECONDS) 

C  CONSTl   =   (43200  •   0.01/1(2.751  -  DENSITY)   •  0.48)1/H»*2 

CONSTl   =  900.0/(12.751   -  DENS  I Ty ) «H»H I 

C  THE   MINIMUM  HATER  EQUIVALENT   WHICH  HILL   ACHIEVE   STABILITY  USINO 

C   THE   ABOVE   DENSITY  FUNCTION   IS  4.7  INCHES 

IFICONSTl  -  0.5)  20tl0.10 

C  THE   MODEL   IS  UNSTABLE  -   INDICATE   THAT    IT    IS   NOT   READY   FOR   USE  NOW. 

C   (IT  MAY  BE   INITIALIZED  AGAIN  BY  AN  OBSERVED  PACK  TEMPERATURE  CARD 

C   AND  STABILITY  HILL  BE  ASCERTAINED  FROM  THE  HATER  EQUIVALENT  AT 

C  THAT  TIME) 

10  OREAOY   =  0 

RETURN 

C  GET  THE   SECOND  CONSTANT 

20  C0NST2  =  1.0  -  CONSTl  -  CONSTl 

C  PERFORM  THE   SIMULATION   IN  TWO  PARTS   (ONE   FOR   EACH   12  HOUR  PERIOD). 

C   -SIMTEMl-  HOLDS  THE  THREE  TEMPERATURES  FROM  THE  PREVIOUS  INTERVAL 

C   THAT   ARE  NEEDED   TO   SIMULATE   SIMTEM2.    THE  NODE   AT   THE   CENTER  OF 

C   THE   PACK.      SIMULATE   THE  FIRST   12   HOURS  NOW 

SIMTEM2  =  (CONSTl  •  ISIMTEMKll  ♦  SIMTEMII3))I  +  (C0NST2  *  SIMTEMl 
112)  ) 

C  THE  AVERAGE   SNOHPACK  TEMPERATURE   IS  THE  AVERAGE  OF  THE  2  NODES 

C   (MIDDLE  AND  GROUND)    IN  BOTH   INTERVALS.     GROUND  TEMPERATURE  IS 

C   CONSTANT,    SO  START   THE   AVERAGE  NOW 

$IMTEM3  =  SIMTEM1I3)   ♦  SIMTEM1I3)   ♦  SIMTEM2 

C  RESET  -SIMTEMl-  TO  THE   TEMPERATURES  OF  THE   INTERVAL  JUST  SIMULATED 

C   FOR   USE   IN  THE   SECOND   12   HOUR    INTERVAL   SIMULATION.      THE  SURFACE 

C   AIR   TEMPERATURE    Is   SPLIT    INTO   A  LOW  AVERAGE    (    IHE4NtHIN)/2    >  AND 

C   A  HIGH  AVERAGE   (    (MEAN+MAX)/2    )   FOR  USE   WITH   THE   TWELVE  HOUR 

c   Intervals,    use  the  low  average  now 

SIMTEMim   =   AMINl    (0.0,  (  (  (TEMPMIN-32.0  l*0.5555555556)*AVETEMC  )/ 
1  2.0) 
SIMTEM1(2)   =  SIMTEM2 

C  SIMULATE  THE  SECOND  12  HOURS  AND  COMPUTE  THE  AVERAGE  SNOHPACK 

C   TEMPERATURE 

SIMTEM2  =  (CONSTl  »  ISIHTEMl(l)  ♦  SIMTEM1I3)))  ♦  (C0NST2  *  SIMTEMl 
1(2)) 

SIMTEM3  =   (SIMTEM3  ♦  SIMTEM2)/4.0 

C  RESET  -SIMTEMl-  USING  THE  HIGH  AVERAGE  FOR  USE  ON  THE  FIRST 

C   INTERVAL  OF  THE  NEXT  DAY 

SIMTEMKl)   =  AMINl   (0.0,  (  (  I  TEMPMAX-32.0  )*0.5555555556)+AVErEHC  1/ 
I  2.0) 
SIMTEH1(2)   =  SIMTEH2 

C  CHECK  TO  SEE   IF  THE  GROUND  TEMPERATURE  SHOULD  BE  RAISED 

IFISIMTEM3  ♦   1.5)  60,40,30 
30   IF(SIMTEM3   ♦  0.51  40,50,50 

40   IF ( SIMTEMl ( 31 .LT. -0.51    SIMTEMl(3)   =  -0.5 

RETURN 
50  SIMTEMl (3)   =  0.0 
60  RETURN 

END 


Subroutine  EVTRAN 


SUBROUTINE  EVTRAN 

C  COMPUTE   THE  EVAPORATION  AND  TRANSPIRATION  DURING  THE  GROHING 

C   SEASON 

C  

C  DICTIONARY 

C  AVABLE  -   THE   FACTOR  FOR   ADJUSTING  THE   EVAPOTRANSP I R AT  I  ON  FOR 

C  AVAILABLE    SOIL  HATER 

C  CANREF   -   THE   FACTOR  FOR  ADJUSTING  THE   EVAPOTRANS P I R AT  ION 

C  FOR  CANOPY  REFLECTIVITY 

C  

COMMON /ONLYCOR/   A VE TEMC , BASTEMF , C ALDEF , COM AX , COVDEN , ORE ADY , ENGB AL , 

1  ENGBALl ,FREEHAT, LAS TUSD,NOAYSN0,ONTR EES, PHASE, PRE WEOV.RECHRG, 

2  SIMTEMl (3 ) ,SIMTEM3,TC0EFF,THRSHL0,VEGTYPE 
INTEGER   OREAOY, PHASE, VEGTYPE 

CDMM0N/WATRBAL/ETFROM,EVAPQTR,GENR0,PRECIP,RADIN,RA0LWN,RADSW»J, 
1  TEMPMAX,TEMPHIN,HATERIN 
ETFROM  =  4.0 

C  GET  THE  ADJUSTMENT  FACTOR  FOR  AVAILABLE  SOIL  WATER.     IF  THE 

C   RECHARGE  REQUIREMENTS  ARE  SATISFIED  (I.E.,   FIELD  CAPACITY  HAS 

C   BEEN  REACHED),   MAXIMIZE   THE   ADJUSTMENT  FACTOR 

IF(RECHRG)  20,20,10 
10  AVABLE   =  1.0 
GO  TO  80 

C  in  FORESTED  AREAS,   USE  THE  MAXIMIZED  FACTOR  UNTIL  THE  RECHARGE 

C   REQUIREMENTS   ARE  HALF  OF   THE  FIELD  CAPACITY.      AT   THAT   POINT,  USE 

G   THE  LINEAR  FUNCTION  Y  =  MX  »  B,  WHERE  B  =  0.0  TO  GET  THE  FACTOR 

20  IFICOVDENI  50,50,30 

30  IF(RECHRG  ♦  2.65)  40,10,10 

40   AVABLE   =   0.377  *   (5.3  ♦  RECHRG) 
GO  TO  8  0 

C  IN  CLEARINGS,   ASSUME   THAT   EVAPOT R ANS P I R AT  I  ON  DECREASES  LINEARLY 

C   FROM  THE   MAXIMIZED  VALUES   AT   FIELD  CAPACITY   TO  ZERO  AT 

C   THREE-FOURTHS  OF  FIELD  CAPACITY 

50   lF(RECHRt.  *   1.325)  70,70,60 

60  AVABLE  =   (0.755  »   (5.3  +  RECHRG))   -  3.0 

GO  TO  80 
70  EVAPOTR  =>  0.0 
RETURN 

C  COMPUTE  THE  FACTOR  FOR  ADJUSTING  THE  EVAPOTRANSPIRATION  FOR  CANOPY 

C   REFLECTIVITY   (PROTECT   AGAINST  DIVISION  BY  ZERO  -  JUST  DEFINE  THE 

C   FACTOR  FOR  CLEARINGS  AS  THE  MINIMUM  VALUE) 


80   IF(CDMAXI  90,90,100 
90   CANREF   =  0.5 
GO   TO  130 

C  COMPARE  THE  COVER  DENSITY  WITH  THE  MAXIMUM  (CHECK  FOR  THINNING) 

100   IFICOVOEN  -    (CDHAX/3.0))  110,110,120 

C  FOR  COVER   DENSITIES   THINNED  TO  ONE-THIRO  OF   THE  MAXIMUM  OR  LESS, 

C   USE  THIS  RELATIONSHIP 

C  CANREF   =  1.0  -   (0.5  -   ( (0.75«CD)/CDMX) ) 

110  CANREF   =  0.5  ♦    ((0.75   »  COVDEN )/ COM AX ) 
GO   TO  130 

C  FROM   MAXIMUM  COVER   DENSITY   DOWN   TO  ONE-THIRO  OF   THAT  VALUE,  USE 

C   THE   FOLLOWING  RELATIONSHIP 

C  CANREF  =  1.0  -   (0.25  -   ( 10. 15/(0. 67«C0MX) )«(C0-(0.33»CDMX) )) ) 

120  CANREF   =   0.75   *   ((0.15/10.67   *  CDMAX ))»( COVDEN  -    (0.33  •  COHAX))) 

C  PERFORM   THE  ADJUSTMENTS 

130  EVAPOTR   =  EVAPOTR   •   AVABLE   »  CANREF 

C  SEE    IF    THE   EVAPOTRANSPIRATION  WILL   DEPLETE   THE   MANTLE  STOARGE 

C   BELOW   THE   WILTING  POINT.      IF   SO,   ALTER   THE  EVAPOTRANSPIRATION 

C   ACCORDINGLY 

IF(RECHRG  -   EVAPOTR  ♦  5.3)  140,150,150 
140    EVAPOTR    =    RECHRG   +  5.3 
RECHRG    =  -5.3 
RETURN 

150   RECHRG    =   RECHRG   -  EVAPOTR 
RETURN 
END 


Subroutine  LINK 


SUBROUTINE   LINK    I C AL AI R ,C ALOR I E , I  RETURN ) 

C  THIS   SUBROUTINE   IS   THE    INTERFACE   BETWEEN   THE  RADIATION  BALANCE 

Q   (SUBROUTINE  RAOBAL)    AND   THE   DIFFUSION  MODEL    (SUBROUTINE  OIFMOD) 

C  

C  DICTIONARY 

C 

C  CALDM  -   THE  CALORIC  LOSS  OF  GAIN  AS  COMPUTED  BY  THE  DIFFUSION 

C  MODEL 

C 

COMMON /ONLYCOR/   AVET EMC , B AST EHF , C AL DE F , COM AX , CO VDEN , DREADY , ENGB AL , 

1  ENGB AL 1, FR EE W AT, L AS TUSD,N0AYSNO,ONTREES, PHASE, PRE WEOV.RECHRG, 

2  SIMTEMl  (3), SI MTEM3,TC0EFF,THRSHL0. VEGTYPE 
INTEGER  OREAOY. PHASE, VEGTYPE 

C OHHON/ WAT RB AL/ ETFROM, EVAPOTR, GENRO,PREC IP, RAD  I N,R AOL WN,R ADS WN, 
1  TEMPMAX.TEMPHIN.WATERIN 

C  SEE   IF  THE  RADIATION  BALANCE   IS  AN  ENERGY  LOSS  OR  GAIN 

IF(CALORIE)  10,10,80 

C  THERE   WAS  A  LOSS.      IF  THIS   IS  STILL  WINTER   (NO  FREE  WATER),  JUST 

C   GO  AHEAD  AND   USE   THE  DIFFUSION  MODEL 

10   IFIFREEWAT)  20,20,50 

C  USE   THE   DIFFUSION   MODEL   TO   SIMULATE   THE   CURRENT   AVERAGE  SNOWPACK 

C  TEMPERATURE 

20   I F (DREADY. NE. 1 )   GO  TO  140 
CALL  OIFMOD 
IFIOREAOYI  40,40,30 

C  NOW  MAKE   ANY  NECESSARY   ADJUSTMENTS   IN  THE  RADIATION  BALANCE  TO 

C   CAUSE   THE  PACK  TEMPERATURE  TO  BE  THE  SAME  AS  -SIMTEM3-.  GET  THE 

C   DIFFERENCE  BETWEEN  THE  CALORIE  DEFICITS  AS  COMPUTED  BY  THE 

C   DIFFERENT  METHODS 

30  CALDM  =  CALDEF  ♦   ISIMTEH3  •  PREHEQV  »  1.27) 

C  ADJUST   THE   LONG  HAVE   PORTION  OF   THE   RADIATION   BALANCE   BY  THE 

C   DIFFERENCE  BETWEEN  THE  CALORIES  DERIVED  FROM  THE  DIFFUSION  MODEL 

C   AND  THE  ENERGY  BALANCE 

CALORIE   =  CALDM 
RAOLWN  =  CALORIE  -  RADSWN 
40    IRETURN  =  0 
RETURN 

C  THE  LOSS  IS  USED  TO  FREEZE  PART  OR  ALL  OF  THE  FREE  WATER,   BUT  IT 

C-   HAY  NOT  CREATE  COLO  CONTENT.      IF   IT  WOULD  CREATE  COLD  CONTENT, 

C   RE-INITIALIZE  THE  DIFFUSION  MODEL  TO  0  AND  ADJUST  THE  ENERGY 

C   BALANCE  ACCORDINGLY 

50  CALL   CALDSS  (CALORIE] 

IFIFREEWAT  -  0.05)  60,60,70 
60  SIMTEMKl)   =  AMINl   (  AVETEMC  ,  0.  0 1 
SIMTEMl(2)   =  0.0 
SIMTEM1(3I    =  0.0 
DREADY  =  1 

C  MAKE  ANY  NECESSARY  ADJUSTMENTS  TO  THE  ENERGY  BALANCE  TO  COMPENSATE 

C   FOR  THE  COLD  CONTENT  THAT  WOULD  HAVE  BEEN  GENERATED  BY  THIS  LOSS 

C   AND  ZERO  THE   COLO  CONTENT 

ENG6AL  =  ENGBAL  *  CALDEF 

RAOLWN   =  RAOLWN  ♦  CALDEF 

FREEWAT   =  0.0 

CALDEF   =  0.0 
70   IRETURN  =  1 

RETURN 

C  THERE   IS  CALORIC   INPUT  TO  THE  PACK.     CHECK  TO  SEE   IF  CONDITIONS 

C   INDICATE   THAT  THE   DIFFUSION  MODEL   SHOULD  BE  TURNED  OFF   AND  THE 

C   ENERGY  BALANCE  USED  FOR  SPRINGTIME   SIMULATION.     CONSIDER  FIRST 

C   ANY  COLD  CONTENT    (INCLUDING   THAT  OF   THE   PREVIOUS   DAY   AND  ANY 

C   CREATED  BY  A  SNOW  EVENT  ON  THIS  DAY).     IF  THERE   IS  COLO  CONTENT, 

C   CHECK   THE   AVERAGE   AIR   TEMPERATURE   AND  THE   SNOwPACK  TEMPERATURE 

C   FROM   THE   PREVIOUS   DAY   FOR   ARBITRARILY  CHOSEN  SPRINGTIME 

C   CONDITIONS  AND  IF  ALL  ARE  NOT  SATISFIED,   GO  AHEAD  AND  USE  THE 

C  DIFFUSION  MODEL 

80   IFICALDEF)  170.170,90 
C  0.889  =   1.27   •  0.7  DEGREES  C   (ARBITRARY  TEMPI 

90   IFIAVETEMC.LE. 0.0. OR. CALDEF. GT.PREWEOV»0. 889)   GO  TO  20 

C  SINCE   SPRINGTIME   CONDITIONS   PREVAIL,   RECOMPUTE   THE   BACK  RADIATION 

C   AND   THE   NET  RADIATION   BALANCE    (REMEMBER,    IF   THERE    IS   SNOW,  THE 

C-   LONGWAVE   IS  ASSUMED  TO  BE  ZERO,   SO  THERE  WOULD  BE  NO  NEED  To  HAKE 

C   ANY  ADJUSTMENTS) 

IF(NOAySNO)  140,140,100 
100  USE  =   (TEMPMIN  -  32.0)   •  0.5555555556 
IFIUSE.GT.O.OI   USE  =  0.0 
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CALSNOW  =   1.17E-7  »    {(USE   ♦   273.  15)  <.  I 

IFIPRECIPI  110,110,120 
110  R60LWN   =    lll.O  -  COVOEN)    •   ((0.  757   •   C4L&1RI   -  C4LSN0HII    <■  (COVDEN 
1   «   (CALAIR  -  CALSNOH) I 
Gd   TO  130 
120   RAOLWN    =   CALAIR   -  CALSNOW 
130  CALORIE   =   RADSHN  »  RAOLWN 

C  RE-INI r lALI ZE   THE  DIFFUSION  MODEL  TO  THESE  CONDITIONS   (BUT   IF  THE 

C   INPUT    IS  MORE   THAN  ENOUGH  TO  WIPE   OUT   THE  CALORIE   DEFICIT,  JUST 

C   LET   IT  BRING  THE   PACK  TO   ISOTHERMAL.      IN  THIS  WAY,   TWO  CONSECU- 

C   TIVE   DAYS   OF    INPUT   ARE   REQUIRED  TO  GENERATE   FREE  WATER) 

140  COMPARE  =  CALORIE  -  CALOEF 
IF(COMPARE)  160,150,150 

C  INITIALIZE   THE   DIFFUSION  MODEL   TO   ISOTHERMAL  CONDITIONS 

150   SIMTEHU  1  )    =  0.0 
SIHTEMl(2)    =  0.0 
SIHTEHl (3)    =  0.0 
SIMTEM3   =  0.0 
OREADY   =  1 
GO  TO  30 

C  REDEFINE  THE  SURFACE  TEMPERATURE  AND  COMPUTE  THE  NEW  AVERAGE  PACK 

(;  TEMPERATURE.      THEN  COMPUTE   THE   MIDDLE   NODE   AS   A  FUNCTION  OF  THAT 

C   AVERAGE,    THE   SURFACE   TFMPERATURE  AND  THE  GROUND  TEMPERATURE 

C   (WHICH    REMAINED  UNCHANGED) 

160    SIMTEMl(l)    =    AMINl    ( 0 . 0 , AVE T EMC ) 

SIMTEM3   =  COMPARE/(PREWEOV   •  1.271 

SIMTEMl(2)   =   (3.0  •  SIMTEM3I   -  SIMTEMl(l)  -  SIMTEM1(3) 
SIHTEM1(3)   =  0.0 
DREAOY  =  1 
GO  TO  30 

C  THERE   IS  INPUT  TO  THE  PACK  AND  THE  PACK   IS   ALREADY   ISOTHERMAL.  IF 

C   THIS  ENERGY  WILL   CREATE   AT  LEAST   0.05   INCH   (ARBITRARY  AMOUNT)  OF 

C   FREE   WATER,   SET   THE   DIFFUSION  MODEL   TO  STANDBY   STATUS  AND   LET  THE 

C   ENERGY  BALANCE  TAKE   ITS  COURSE 

170   IF(FREEWAT   ♦    ( CA LOR  I E/ 203 . 2 )    -  0.05  )    150,  180,  180 
180  DREADY  =  0 
IRETURN   =  0 
RETURN 
END 


1   .64,    .63,    .62f   .61,  .60/ 
DATA  REFMELT/.72t   .65,   .60,   .58,   .56,   .5*,    .52,    .50.   .48,  .4 
1    .44,    .43,    .42,    .41,  .40/ 
PASTINT   =  NDAYSNO 
IF(NDAYSNO)  80,80,10 

C  USE   THE    SAME   FUNCTION  AS   LAST  TIME 

10   IF(LASTUSD)  20,20,50 

C  ACCUMULATION  PHASE  -  AFTER  15  DAYS,  USE  THE  MELT  FUNCTION 

C   STARTING  AT   THE   FOURTH  DAY 

20   IF(PASTINT  -   15)  30,30,40 
30  PACKREF   =  REFACUM(PASTINT) 
RETURN 

40  PASTINT   =  PASTINT  -  II 

C  — MELT  FUNCTION  -  AFTER   15  DAYS.  USE  A  CONSTANT  40  PERCENT 

50   IF(PASTINT  -   15)  70,70.60 
60  PASTINT   =  15 
70  PACKREF   =  REFHELT(PASTINT) 
RETURN 

C— - — THERE   IS  NEW  SNOW  -  DETERMINE   IF  THE  FUNCTION  IS  TO  BE  RE- 
C— - —  INITIALIZED 

60   IFITEMPHAX  -   THRSHLD)  90.90.10 

C  IT   IS,   SO  SEE  WHICH  FUNCTION   IS  TO  BE  USED 

90   IFICALDEFI  110,110,100 
100  PACKREF   =  0.91 
LASTUSD   =  0 
RETURN 

C  THE  PACK   IS  ISOTHERMAL.   BUT  IF  THE  ENERGY  BALANCE  FROM  THE 

C   PREVIOUS  INTERVAL  WAS  NEGATIVE,  USE  THE  ACCUMULATION  PHASE 

C   FUNCTION  ANYWAY 

110  IF(ENGBALl)  100.120.120 
120  PACKREF  =  0.81 
LASTUSD   =  1 
RETURN 
END 


Subroutine  MIXTURE 


SUBROUTINE  MIXTURE 

C  THIS   SUBROUTINE   CONTROLS   THE   COMPUTATIONS  FOR   A  PRECIPITATION 

C   EVENT  THAT   IS  A  MIXTURE  OF  SNOW  AND  RAIN 

C  

COMMON/ ONLYCOR/  AVE T EMC , BA STE MF , CAL DE F , COM AX , COVDEN , ORE ADY , E NGBAL , 

1  ENGBALl, FREE  WAT, LASTUSD, NO A YSNO,ONTREES,PHASE,PREWEOV,BECHRG, 

2  SIMTEMl (3) , SI MTEM3, TCOEFF , THRSHLD, VEGTYPE 
INTEGER  DREADY, PHASE, VEGTYPE 

COMMON/WATRBAL/ETFR0M,EVAPOTR,GENRO,PRECIP,RAOIN.RADLWN,RADSWN, 
1  TEHPMAX,TEMPMIN,WATERIN 
C  DICTIONARY 


C  AMTRAIN  -  THE  AMOUNT  OF  PRECIPITATION  OCCURRING  AS  RAIN 

C  TFORAIN  -   THE   TEMPERATURE  FOR  COMPUTING  THE  DEPLETION  OF  THE  TOTAL 

C  CALORIE  DEFICIT   CAUSED  BY   THE   RAIN   (DEGREES  CI 

C  TFORSNU  -   THE   TEMPERATURE   FOR  COMPUTING  THE   CONTRIBUTION  OF  THE 

C  SNOW  TO  THE  TOTAL  CALORIE  DEFICIT   (DEGREES  CI 

C  

C  COMPUTE   THE  AMOUNT  OF  PRECIPITATION  OCCURRING  AS  RAIN 

C   AMOUNT   RAIN  =   P   •   ( B/ A ) ,  WHERE 

Q   p   -   PRECIPITATION   IN  INCHES 

C   B  =  DAILY  MAXIMUM  TEMPERATURE  -  BASE  TEMPERATURE   (DEGREES  F) 

C   A  =  DAILY  MAXIMUM  TEMPERATURE  -  MINIMUM  TEMPERATURE   (DEGREES  Fl 

B  =   TEMPMAX  -  BASTEMF 

A   =    TEMPMAX   -  TEMPMIN 

AMTRAIN   =  PRECIP   «  (B/A) 

C  NOW  COMPUTE  THE  AVERAGE   TEMPERATURES   (DEGREES  C)  WHICH  PRODUCE 

C-  SNOW  AND  RAIN 

TFORSNO  =    (TEMPMIN  ♦   BASTEMF   -  54.0)    *  0.2777777778 

TFORAIN  =   (TEMPMAX  +  BASTEMF  -  54.0)   »  0.2777777778 
C  COMPUTE   THE   EFFECT   OF   THE   SNOW  ON  THE  SNOWPACK 

CALL   SNOWED   I TFORSNO, PRECI P-AMTRAIN) 

C  COMPUTE    THE   EFFECT  OF   THAT   PORTION  OF   THE   PRECIPITATION  OCCURRING 

C   AS  RAIN  ON  THE  SNOWPACK 

CALL   RAINED   I TFOR A  I N, AMTRA I N ) 

RETURN 

END 


Function  PACKREF 


FUNCTION  PACKREF  (DUMMY) 
—GET   THE   REFLECTIVITY  OF  THE  SNOWPACK 

COMMON /ONLYCOR/  AVETEMC , BA STEMF . CALDEF , COMAX . COVDEN . ORE ADY , ENGBAL . 

1  ENGBALl .FREE  WAT, LASTUSD, NOAYSNO.ONTREES, PHASE, PREWEQV.RECHRG, 

2  SIMTEMl (3) ,SIMTEM3,TC0EFF,THRSHL0,  VEGTYPE 
INTEGER  DREADY, PHASE. VEGTYPE 

COMMON / WATRBAL/E TF ROM, EVAPOTR.GENRO, PRECIP, RAO  I N,R AOL WN,R ADS WN, 
1   TEMPMAX, TEMPMIN, WATERIN 


C  DICTIONARY 


PASTINT  -  A  VARIABLE  SET  EQUAL   TO  -NDAYSNO-  AND  ALTERED  AS  NEEDED 
TO  CHOOSE  THE  PROPER  REFLECTIVITY  FUNCTION 

REFACUH  -  A  REFLECTIVITY  FUNCTION  FOR   THE   SNOWPACK   DURING  THE 
ACCUMULATION  PHASE  OF   THE  SNOWPACK 

REFHELT  -  A  REFLECTIVITY  FUNCTION  FOR   THE  SNOWPACK   DURING  THE 
MELT  PHASE  OF  THE  SNOWPACK 

DIMENSION  REFACUM(15) ,REFMELT(15I 
INTEGER  PASTINT 

DATA  REFACUH/. 80.    .77,    .75,    .72.    .70.    .59,    .68,    .67.    .66,  .65, 


Subroutine  RADBAL 


SUBROUTINE  RADBAL 

THIS  SUBROUTINE  COMPUTES  THE  RADIATION  BALANCE  AND  TRANSFERS 

CONTROL  TO  THE  DIFFUSION  MODEL   IF   IT  IS  NEEDED 
COMMON /ONLYCOR/  AVETEMC, BASTEMF, CALDEF. COM AX. COVDEN. ORE ADY, ENGBAL. 

ENGBALl , FREE  WAT, LASTUSD, NOAYSNO.ONTREES. PHASE, PRE WEQV,RECHRG, 

SIMTEMl (3) ,SIMTEH3,TC0EFF, THRSHLD. VEGTYPE 
INTEGER  DREAOY, PHASE. VEGTYPE 

COMMON / WATRBAL/E TF ROM, EVAPOTR.GENRO. PRECI P. RAD  I N,R AOL WN,R ADS WN, 
TEMPMAX, TEMPMIN, WATER  IN 


C  DICTIONARY 


C  CALAIR  -  POTENTIAL  LONGWAVE  CALORIC  INPUT  AT  AIR  TEMPERATURE 

C  CALORIE  -  CALORIES  OF  HEAT  ABSORBED  OR  RELEASED  BY  THE  SNOWPACK 
C  FROM  THE  NET  RADIATION  BALANCE 

C  CALSNOW  -   POTENTIAL   LONGWAVE  CALORIC  LOSS    AT   SNOW  TEMPERATURE 

C  SNOCAN  -  THE  LONGWAVE  RADIATION  BALANCE  BETWEEN  THE  SNOW  AND  THE 
C  CANOPY 

C  SNOSKY  -  THE  LONGWAVE  RADIATION  BALANCE  BETWEEN  THE  SNOW  AND  THE 
C  SKY 

C  COMPUTE   THE  CALORIC    INPUT   FROM  NET   SHORT  WAVE  RADIATION  AS  A 

c   Function  of  the  snowpack  reflectivity 

RADSWN  =  RADIN  *    (1.0  -   PACKREF    (0.01)   •  TCOEFF 

C  IF  THE  PACK   IS  ACCUMULATING.   BUT   IS  NOT  DEEP  ENOUGH  FOR  STABILITY 

Q   IN  THE  DIFFUSION  MODEL,  USE  THE  FOLLOWING  SIMPLIFIED  METHOD  FOR 

C   DERIVING  THE   RADIATION  BALANCE 

IF(PHASE)  60,10,100 
10  IFIPrEWEQV  -  4.7)  20.50.50 

C  USE  ONLY  THE  SHORTWAVE   INPUT   (THIS  IMPLIES  THAT  THE  ONLY  COLD 

C   CONTENT  GENERATED  IN  THE  ACCUMULATING  PACK   IS  THAT  OF  NEW  SNOW) 

20  CALORIE   =  RADSWN 

RAOLWN   =  0.0 

CALL  CALIN  (CALORIE) 

C  IF  THE  MEAN  TEMPERATURE  WAS  LESS  THAN  OR  EQUAL  TO  0  C,   00  NOT 

C   ALLOW  ANY  MELT  OR  FREE  WATER 

IFIAVETEHC)  30,30,40 
30  PREWEQV  =  PREWEQV  ♦  WATERIN 

RAOLWN   =  -ENGBAL 

ENGBAL   =  0.0 

WATERIN  =  0.0 

FREEWAT  =  0.0 

CALDEF  =  0.0 
40  RETURN 

C  THE  PACK  HAS  JUST  REACHED  A  SUFFICIENT  DEPTH.     INITIALIZE  THE 

C   DIFFUSION  MODEL,   BUT  RETAIN  PSEUDO-CONTROL  UNTIL  THE  DIFFUSION 

C   MODEL   IS   WELL   ALONG   INTO  STABLE  CONTROL 

50  PHASE  =  -1 
DREADY   =  1 

C  START  THE  PACK  AT  -3  C   (CAL  DEF  =  PACK  TEMP  *  PREWEQV  *   1.27  1 

C  3.81  =  3  »  1.27 

CALOEF  =  3.81  »  PREWEQV 

RAOLWN   =  -   CALOEF   -  RADSWN  -  ENGBAL 

ENGBAL   =  -  CALDEF 

SIMTEMKl)    =  AHINl    (  AVETEMC  .  0  .  0  I 

SIMTEMl ( 2)   =  -3.0 

SIMTEM1(3)    =  -1.5 

FREEWAT   =  0.0 

RETURN 

C  THE   DIFFUSION  MODEL   HAS   BEEN   INITIALIZED  PREVIOUSLY.      IF   IT  IS 

C   STILL   STABLE  AND   IF  THE  PACK  IS  DEEP  ENOUGH  TO  INSURE  CONTINUED 

C   STABILITY  UNTIL  MELT,   RELINQUISH  CONTROL   COMPLETELY   TO  THE 

Q   NORMAL   METHOD  OF   COMPUTING  THE   RADIATION  BALANCE.    INTERFACED  WITH 

C  THE   DIFFUSION  MODEL 

50   IF(DREADY)  100,70,80 

70  PHASE   =  0 
GO  TO  10 

80  CALORIE  =  RADSWN 
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CALAIR  =  0.0 
RADLWN  =  0.0 

IFIPREwEav  -  5.01  160,160,90 
90  PHASE  =  1 

C  USE   THE   NORMAL   METHOD  OF  COMPUTING   THE   RADIATION   BALANCE.      IF  ANY 

C   OF   THE   PRECIP  WAS   SNOW,   THE   NET   LONG  WAVE   RADIATION   BALANCE  IS 

C  ASSUMED  TO  BE  ZERO 

100  IF(NOAVSNO)  110,110,120 
110  RADLWN   =  0.0 
CALAIR   =  0.0 
GO  TO  150 

C  TO  COMPUTE   THE  LONG  WAVE  RADIATION  COMPONENTS,   CONVERT   THE  AIR 

c   AND  Snow  temperatures  to  potential  calories  by  the  stefan  - 

c  blotzmann  function,  calories  ■=  s  »  it       '.I,  where 

C  S  =   1.17E-7  CAL/I  (CH*»2)  (DEGREES  KELV  IN  )••'>) ,  AND 

C   T  =  ABSOLUTE   TEMPERATURE   (DEGREES  KELVIN) 

120  CALAIR   =   1.17E-7   »    ((AVETEMC  ♦   273.16)   •»  4) 
USE  =  AVETEMC 

C  IF  THE    SNOWPACK   IS   ISOTHERMAL,   USE   THE   MINIMUM  TEMPERATURE  FOR 

C   COMPUTING   THE   BACK  RADIATION 

IF (CALDEF.EQ.0.0)   USE  =   (TEMPMIN  -  32.0)   »  0.5555555556 

C  UNDER  NO  CIRCUMSTANCES  MAY  THE  TEMPERATURE  FOR  COMPUTING  THE  BACK 

C   RADIATION  BE  GREATER  THAN  ZERO 

IF(USE.GT.O.O)   USE  =  0.0 

CALSNOW  -  1.17E-7  *   ((USE  »  273.16)   *»  4  1 

C  COMPUTE   THE   LONG  WAVE  RADIATION  COMPONENTS  AS  A   FUNCTION  OF  THE 

C   FIRST,   DETERMINE  WHETHER  THE  SKIES  ARE  CLEAR  OR  CLOUDY 

IF(PRECIP)  130,130,140 

C  WITH  CLEAR   SKIES,    THE   DOWNWARD  LONGWAVE  RADIATION  COEFFICIENT  IS 

C  .757    (RUNOFF   FROM  SNOWMELT,    EM  1 1 10-2- 1406 ,    US   ARMY  CORPS  OF 

C   ENGINEERS,   1960,   PAGE  7) 

130  SNOSKY   =   (1.0  -  COVDEN)   •   ((0.757  «  CALAIR)   -  CALSNOW) 

C  THE   DOWNWARD  LONGWAVE   RADIATION  COEFFICIENT    IS    l.O   BENEATH  THE 

C   FOREST  CANOPY   (OR  BENEATH  CLOUDY  SKIES) 

SNOCAN  =  COVDEN  »    (CALAIR  -  CALSNOW) 

RADLWN   =   SNOCAN  *  SNOSKY 

GO  TO  150 

C  WITH  CLOUDY   SKIES,   WHEN   THE   DOWNWARD  LONGWAVE  RADIATION  COEFFI- 

C   CIENT  IS  1.0  INSTEAD  OF  .757,   THE  ABOVE  THREE  EQUATIONS  MAY  BE 

C   REDUCED  ALGEBRAICALLY  TO  THE  FOLLOWING  SINGLE  EQUATION 

140  RADLWN  =  CALAIR  -  CALSNOW 

C  COMPUTE   THE   CALORIC   INPUT  OR   LOSS  FROM  THE  NET   EFFECT  OF  SHORT 

C   WAVE   AND  LONG  WAVE  RADIATION 

150  CALORIE   =   RADSWN  ♦  RADLWN 

C  THE  SNOWPACK  TEMPERATURE  DIFFUSION  MODEL   (LEAF,    1970,   STUDY  PLAN 

C   FS-RM-1602,   NO.   224.     ROCKY  MOUNTAIN  FOREST   ANO  RANGE   EXP   STA)  IS 

C   INCORPORATED   TO  CONTROL  THE  SNOWPACK  TEMPERATURE  AND  COLD  CONTENT 

C   DURING  NON-ISOTHERMAL  CONDITIONS.     SEE  NOW   IF  THE  DIFFUSION  MODEL 

C   MAY  BE   USED   (DREADY  MAY  NOT  BE   -1   AND   PASS   THROUGH  LINK   SINCE  IT 

C   IS  NOT  DESIGNED   TO  WORK   WITH   IT.     THE  -1    IS   USED  TO   INDICATE  THAI 

C   THE  RADIATION  ROUTINES   ARE  TO  BE  USED  EXCLUSIVELY).     IF   IT  MAY  BE 

C   USED,   PASS  THROUGH  THE  LINKING  ROUTINE  WHICH   INTERFACES  THE 

C   DIFFUSION   MODEL   AND  THE   RADIATION  ROUTINES 

IF(DREADY)  170,160,160 

160  CALL   LINK    (CALAIR, CALORIE, IRETURN) 
IF(IRETURN)  170,170,190 

170  IF(CALORIE)  180,190,200 

180  CALL  CALOSS  (CALORIE) 

190  RETURN 

200  CALL  CALIN  (CALORIE) 
RETURN 
END 


Subroutine  RAINED 


SUBROUTINE   RAINED  (TFORAIN,AMTRAIN) 

C  THIS  SUBROUTINE  COMPUTES  THE  EFFECT  OF  RAIN  ON  SNOW 

COHMON/ONLYCOR/   AVE TEMC , B AS TEMF , CALDEF , COM AX , COVDEN , DREADY , ENGBAL , 

1  ENGBAL 1 ,FREE WAT, L«STuSD,NDAYSNO,0NTREES,PHASE,PREW6QV,RECHRG, 

2  SIMTEMl 13) , SIMTEM3, TCOEFF ,THRSHLD,VEGTYPE 
INTEGER  DREADY, PHASE .VEGTYPE 

COMMON/ WATRB AL / E TFROM, EV AP0TR,GENRO, PRECIP, RAD  I N,R AOL WN,  RADSWN, 
I   TEMPMAX, TEMPMIN, WATERIN 


C  ALL   OF   THE   RAIN   IS   ADDED  TO  THE   FREE   WATER   AND  CONTRIBUTES  CALORIC 

Q  INPUT    TO    THE  PACK 

50   FREEWAI    =   FREEWAT    t  AMTRAIN 

CALL  CALIN    (TFOBAIN   »   AMTRAIN  *  2.54) 

RETURN 

END 


Subroutine  SNOWED 


SUBROUTINE    SNOWED   ( TFORS NO , AMTSNOW ) 

C  THIS   SUBROUTINE  COMPUTES  THE  EFFECTS  OF  A  SNOW  EVENT  ON  THE 

C  SNOWPACK 

COMMON /MASTER /DATE  13) , TMXMSTR , TMNMSTR, PPTHSTR, PPTONOW , OBS HYDR , 
1  POTRAO,MSTREOF,IYR 
INTEGER  DATE 

COMMON /ONLYCOR/  AVE TEMC , BASTE MF , CALDEF , CDM AX , COVDEN , ORE ADY ,  ENGBAL , 

1  ENGBALl ,FRFEWaT,LASTUSO,NOAYSNO,ONTREES,PHASE, PREWEQV,RECHRG, 

2  SIMTEMl 13) ,SIHTEH3,TC0EFF,THRSHLD,VEGTYPE 
INTEGER  DREADY, PHASE, vEGTyPE 

COMMON / WATRB AL/E TFROM, EVAP0TR,GENR0,PREC I P, RAD  in, RADLWN, RADSWN, 
1    TEMPMAX, TEMPMIN. WATERIN 
REAL  INTRCPT 

C  

C  DICTIONARY 

C  

C  AMTSNOW  -  THE   AMOUNT  OF  PRECIPITATION  OCCURRING  AS  SNOW  I  INCHES) 

C  CALSNOW  -   THE   CONTRIBUTION  OF   THIS   SNOW   TO  THE   TOTAL  CALORIE 

C  DEFICIT  ICALORIES) 

C  INTRCPT  -   THE   AMOUNT   OF   SNOW   INTERCEPTED  DURING  THIS  PRECIP  EVENT 

C  TFORSNO  -  THE  TEMPERATURE  FOR  COMPUTING  THE  CONTRIBUTION  OF  THIS 

C  SNOW   TO  THE   TOTAL   CALORIE   DEFICIT    I  DEGREES  C) 

C  

C  DO  NOT   ALLOW  ANY   INTERCEPTION   IN  JULY   AND  AUGUST 

IF (DATE  I  I  I .EQ.7.0R.0ATE I  1) .EQ.8 )   GO   TO  10 

C  DETERMINE   THE   AMOUNT  OF    INTERCEPTED   SNOW  AS   A  FUNCTION  OF  COVER 

C   COMPOSITION   AND  COVER  DENSITY 

IFICDMAX)  10,10,20 
10   INTRCPT   =  0.0 

GO  TO  80 
20  IFIVEGTYPE  -  1)  30,30,40 

C  LODGEPOLE  PINE 

30  PERCENT  =  0.10 
GREATST  =  0.20 
GO  TO  50 

C  SPRUCE  FIR 

40  PERCENT  =  0.15 
GREATST   I  0.30 
50   INTRCPT  =  AMTSNOW  «  PERCENT  •   I COVDEN/COHAX ) 

IFIONTREES  +  INTRCPT  -  GREATST)  70,70,60 
60   INTRCPT  =  GREATST  -  ONTREES 
70  ONIREES  =  ONTREES  +  INTRCPT 
80  NDAYSNO   =  0 

C  ADO   THIS   AMOUNT  OF   PRECIPITATION  TO  THE  PREDICTED  WATER  EQUIVALENT 

PREHEQV   =   PREWEOV   t   AMTSNOW  -  INTRCPT 

C  THE   SNOW  FALLING   WHEN   THE   TEMPERATURE    IS   BETWEEN   35  AND   32  DEGREES 

C   DOES   NOT   ALTER   THE   CALORIC  DEFICIT 

IFITFORSNO.GE.0.0)  RETURN 

C  COMPUTE    THE  CALORIE   DEFICIT   FOR  THIS   SNOW   BY  THE  EQUATION 

C   CALORIE   DEFICIT   -   SII)«OELTA   T»P,  WHERE 

C   S(I)   =   SPECIFIC  HEAT  OF   ICE   (.5  CAL/CM/DEGREES  C). 

C   DELTA  T  =  CHANGE   IN  TEMPERATURE  WITH  RESPECT  TO  FREEZING  (0.0 

C   DEGREES  CENTIGRADE),  AND 

C  p   -  PRECIPITATION   IN  CM   (CONVERSION  FACTOR  =   2.54  CM/ IN  )  . 

C   THEREFORE,   CALORIE  DEFICIT  =  0.5  »  ITfORSNO)   »   (AMTSNOW  •  2.54) 

CALL  CALOSS  (TFORSNO  •   (AMTSNOW  -  INTRCPT)   •  1.27) 

RETURN 
END 


Subroutine  SNOWVAP 


C  DICTIONARY 

C  AMTRAIN  -  THE   AMOUNT   OF   PRECIPITATION  OCCURRING  AS  RAIN  (INCHES) 

C  CALRAIN  -  THE  DEPLETION  OF  THE   TOTAL   CALORIE  DEFICIT  BY   THIS  RAIN 

C  (CALORIES) 

C  TFORAIN  -   THE   TEMPERATURE   FOR  COMPUTING  THE   DEPLETION  OF   THE  TOTAL 

C  CALORIE   DEFICIT   CAUSED  BY   THIS   RAIN   (DEGREES  C) 

C  

C  ADD  THIS  AMOUNT  OF   PRECIPITATION   TO  THE   PREDICTED  WATER  EQUIVALENT 

PREWEQV  =  PREWEOV  ♦  AMTRAIN 
C  SEE   IF  THERE   IS  A  CALORIE  DEFICIT   IN  THE  PACK 

IF(CALOEF)  50,50,10 

C  COMPUTE   THE   AMOUNT  OF   RAIN  AT   THIS   TEMPERATURE   THAT   IS  NEEDED  TO 

C   WIPE   OUT  THE   DEFICIT  AND  COMPARE   IT  WITH  THE  ACTUAL  AMOUNT 

10  CALRAIN  =    (80.0   ♦   TFORAIN)   *  2.54 

AMTNEEO  =  CALOEF/CALRAIN 

COMPARE  =  AMTRAIN  -  AMTNEED 

IF(COMPARE)  30,20,40 

C  THERE  WAS  JUST  ENOUGH  TO  WIPE  OUT  THE  DEFICIT 

20  CALDEF  =  0.0 

ENGBAL   =   ENGBAL  *  CALRAIN 

RETURN 

C  THERE  WAS  NOT  ENOUGH  TO  WIPE   IT  OUT  COMPLETELY.     JUST  DEPLETE 

C  THE  DEFICIT 

30  CALDEF  =  CALDEF  -   (CALRAIN  »  AMTRAIN) 

ENGBAL  =  ENGBAL  »   (CALRAIN  •  AMTRAIN) 

RETURN 

C  THERE   WAS  MORE   THAN  ENOUGH  TO  WIPE   OUT   THE  DEFICIT.     THE  AMOUNT 

C   OF  RAIN  NOT  FROZEN  IS  FREE  WATER 

40  FREEWAT  =  COMPARE 

CALL  CALIN   (TFORAIN  •  COMPARE  •  2.54) 

RETURN 


SUBROUTINE  SNOWVAP 
—COMPUTE   THE  EVAPORATION  FROM  THE  SURFACE  OF  THE  SNOWPACK  AS  A 
—   FUNCTION  OF   The   COVER   DENSITY   AND  REDUCE   THE  PACK  ACCORDINGLY 

COMMON/ONLYCOR/   AVETEMC , BAST EMF , CAL OEF , CDM AX , CO VDEN , DREADY , E NGBAL , 

1  ENGBAL 1, FREE WAT, LA STUSD, NDAYSNO, ONTREES, PHASE, PRE  WE QV,RECHRG, 

2  SIMTEMl (3) , SI MTEM3 .TCOEFF , THRSHLD , V EGT YPE 
INTEGER  DREADY, PHASE, VEGTYPE 

COMMON/ W AT RR AL/ ET F ROM, EVAPOTR.GENRO, PRECIP, RAD  IN, RADLWN, RADSWN, 
I   TEMPMAX, TEMPMIN, WATERIN 
ETFROM    =  2.0 

EVAPOTR  =   11.0  -  COVDEN)   •  EVAPOTR 
PREWEQV  =  PREWEQV  -  EVAPOTR 
RETURN 
END 


Subroutine  ETCODE 


Subroutine  etcode  i etfrom,compsi2) 

C  keep    track  of   WHICH   SOURCES  WERE  USED  FOR  E V APOT R ANS P I R AT  I  ON 

IFIETFROH  -   2.0)  10,20,30 

C  EVAPORATION  FROM   INTERCEPTED  SNOW  ON  CANOPY    lETFROM   =  1) 

10  1FICOMPS12.NE.I.O.AND.COMPS12.NE.3.0.AND.COMPS12.NE.5.0.AND. 
1  C0MPS12.NE.7.0)   C0MPS12  =  C0MPSI2  ♦  ETFROM 
RETURN 

c  Evaporation  from  snowpack  (Etfrom  =  2) 

20   IF ICOMPS12.NE.2.0.AND.COMPS12.NE.3.0.AND.COMPS12.LT.6.0)  C0MPS12 
1  C0MPSI2  ♦  ETFROM 
RETURN 

C  EVAPOTRANSPIRATION  FROM  MELT,   PRECIP  OR   STORAGE   (ETFROM  =  4) 

30  IF(C0MPS12.LT.4.0)   C0MPS12  =  C0MPS12  *  ETFROM 
RETURN 
END 
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Subroutine  GENDATA 


SUBROUTINE   GEND4TA  (N) 

C  GENERATE   THE  DAT*  FOR  THIS  SUBSTATION 

C  

C  DICTIONARY 

C 

C  DO  -  THE  EXACT  POINT  IN  THE  DEGREE-DAY  TABLE  UHICH   IS  TO  BE  USED 

C  IN  THE  COMPUTATION  OF  THE    INCOMING  RADIATION 

C  DOFACT  -  THE  TABLE   OF   ADJUSTMENTS  FACTORS  FOR  COMPUTING  THE 

C  INCOMING  RADIATION 

C  DDI  -  A  REAL.   TRUNCATED  VALUE  OF  -DO-  USED   IN  INTERPOLATION 

C  ET  -   THE   ADJUSTED  POTENTIAL   EV APOTR ANSP IR AT  I  ON    IMAINTAINED  FOR 

C  THE  REDEFINITION  OF  -EVAPOTR-  BY  THE  ALTERNATIVES 

C  RADHORZ  -   THE    INCOMING  RADIATION  COMPUTED  FROM  THE   POTENTIAL  AT 

C  A  HORIZONTAL  SURFACE 

C 

COMMON  AIRTEMC(25,6) 

COMMON  CALDEF(25I ,CDHAX(25) ,COV0EN(25) 
COMMON  DREADY125) 

COMMON  ENGBAL(25) ,ET,ETDAILYI25,12) 

COMMON  FREEWATI25) 

COMMON  IS0THRMI25,20I 

COMMON  LASTUSD (251 .LEVELl ,LEVEL2 

COMMON  MMDO 

COMMON  NDAYSN0I25I ,N0IVSBL,NSU8.NVEARS 
COMMON  QnTREES(25) 

COMMON  PEAKPPTI25.20) , PEAKHE 1 25 . 20 ) . PHASE  I  25 ) .POTENT  I  24) . 
I  PREHE0V(25) 
COMMON  RECHRG(25I 

COMMON  SIMTEHl  (25.31  ,SUBI0(25.6  )  .SLPASP(25.2'.) 

COMMON  TC0EFF(25I .THRSHLDI25I .TOPLOTI 111 

COMMON  VEGTYPE(25I 

COMMON  WEIGHTI25I.HSHEDID(6) 

COMMON    YEARS(20I  .YYMHDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDID 

INTEGER   YEARS, YYMMOO 

COMMON /MASTER /DATE (3  I , TMXMSTR , TMNMSTR , PPTHSTR. PPTONOW , OBSHYDR . 
I   POTRAD.MSTREOF, lYR 
INTEGER  DATE 
COMMON /R AD/ FR AC  TON. I  SUB 

COMMON/WATRBAL/ETFROM.EVAPOTR.GENRO.PRECIP.RAOIN.RADLWN.RADSHN, 
I  TEMPMAX.TEMPHIN.HATERIN 
EQUIVALENCE   I  DATE ( 1 ) , MONTH  I 
DIMENSION  DDFACT(26) 

DATA   DOFACT/   .20.    .35,    .'.S,    .51,    .56,    .59,    .62.    .6*.    .655,  .67. 

1  .682.    .69.    .70,    .71.    .715.    .72.    .722.    .72'..    .726,    .728,  .73, 

2  .734,   .738,   .742.   .746,  .75/ 

C  ADJUST   THE  TEMPERATURES 

IFIDATE(3)  -  671  10,20,20 
10   I   =  5 

GO   TO  30 
20  I   =  3 

30   TEMPMlN  =   AIRTEMC(N,n    ♦   (TMNMSTR  »   A I RTEHC ( N. I ♦! ) ) 
TEMPMSX   =  AIRTEMCIN.ll    ♦    (TMXMSTR  •   A I RTEMC ( N. 2 1 ) 
IF ( TEHPMAX.GE.TEMPMINI   GO  TO  40 
XXX  =  TEMPMAX 
TEMPMAX    =  TEMPMIN 
TEMPMIN    =  XXX 

C  COMPUTE   THE   INCOMING  RADIATION  AT   THE   BASE   STATION   FROM  THE 

C   POTENTIAL   BY  THE  DEGREE-DAY  METHOD 

40  GO  TO   150.50.50.50.60.70.70,70.60,50.50.501 .MONTH 

C  OCTOBER  -   APRIL.   DEGREE   DAYS  =   .44   •   TEMPMAX  -   15.9   (+1.0  FOR 

C   SUBSCRIPTING! 

50   DO  =    (0.44   •   TEMPMAXI   -  14.9 
GO   TO  100 

C  MAY  AND  SEPTEMBER.  DEGREE  DAYS  =  .53  »  TEMPMAX  -  19.5   1+1.0  FOR 

C   SUBSCRIPTING) 

60  00  =    (0.53  •  TEMPMAX)   -  18.5 
GO  TO  100 

C  JUNE.   JULY  AND  AUGUST.   DEGREE   DAYS   =    .63   •  TEMPMAX  -   24.1  (+1.0 

C   FOR   SUBSCRIPTING).   EXCEPT  ON  DAYS  WITH  PRECIP.     DURING  THESE 

C   MONTHS.  USE  A  CONSTANT   44  PERCENT  ON  PRECIP  DAYS 

70  IFIPPTmSTr)  90.90.80 
80  RADHORZ  =  POTRAD  •  0.44 

GO  TO  150 
90  DO  =    (0.63   •  TEMPMAX)   -  23.1 

C  WATCH  FOR   THE   BOUNDARY  VALUES,    0.    AND  25.      (WITH   THE   1.0  ADDED 

C   ABOVE,    THE   SUBSCRIPTS   FOR   THE  TABULAR  VALUES  VARY   FROM   I   TO  261 

100   IF(DO  -    1.0)  110,110,120 

C  USE  THE  FIRST   TABLE  VALUE    (NO  INTERPOLATION  IS  NECESSARY) 

no  RADHORZ   =  POTRAD  •  ODFACT(l) 

GO  TO  150 
120  IF(DD  -  26.0)  140,130.130 

C  USE   THE   LAST   TABLE   VALUE    (NO   INTERPOLATION   IS  NECESSARY) 

130  RADHORZ   =  PCTRAD  *  D0FACTI26) 
GO   TO  150 

C  THE   SUBSCRIPT   IS   IN  THE   PROPER   RANGE.     OBTAIN  THE  INTERPOLATION 

C   FRACTION  AND   SUBSCRIPTS  THROUGH  TRUNCATION  OF  -00- 

140  Jl   -  DD 

DDI  =  Jl 

J  =  Jl   +  I 

C  THE   TERM   (DD-0D1)/1.0   IS   THE   INTERPOLATION  FRACTION 

RADHORZ  =  POTRAD  »  (ODFACT(Jl)  +  ((DDFACT(J)  -  DDFACT(Jll)  •  ( DD 
1   DDI ) ) ) 

C  ADJUST  THE   POTENTIAL  E V APOTR ANS P I R A T I  ON  AS  COMPUTED  BY  THE  HAMON 

C   METHOD  FOR   AVAILABLE   RADIATION   AS   A   PERCENT   OF  POTENTIAL 

150  EvAPOTR  =   ETDAILYIN. MONTH)    •   ( RADHORZ / POTRAD ) 
ET  =  EVAPOTR 

C  ADJUST   THE   RADIATION  AT   THE   BASE  STATION  FOR   SLOPE   AND  ASPECT 


RADIN  =  RADHORZ  »  ( SLPAS P ( N , I  SUB )  +  ( ( SLP ASP ( N, I SUB+1 1  -  SLPASPIN, 
1    ISUal )    *   FRACTONI ) 

C  ADJUST  THE   PRECIP  TO  ENSURE  REACHING  THE  PEAK  WATER  EQUIVALENT 

IF(PPTMSTR)  160.160,170 
160  PRECIP  =  0.0 
RETURN 

170   IF( PEAKPPT(N, lYR I   -   PPTONDHI  190,190,180 

180  PRECIP  =  PPTMSTR  *   ( (PEAKWEIN.IYR)   -  PREhEOV (N ) I / ( PEAKPPI ( N, lYR ) 

1   -   PPTONOHI I 
RETURN 
190  PRECIP   =  PPTMSTR 
RETURN 
END 


Subroutine  PLOTTER 


SUBROUTINE  PLOTTER 

C  PLOT  THE   INFORMATION.     THE  NORMAL  SCALE  IS  20  PRINT  POSITIONS  =  I 

C   INCH,   BUT  SEVERAL  OF  THE  PLOTS  HAVE  ADDITIONAL  SCALE  FACTORS  AS 

C   EXPLAINED  BELOW  TO  ENHANCE  THEIR  VISIBLE  REPRESENTATIONS 

C  

C  DICTIONARY 

C 

C  BOUNDL  -  THE  LOWER  BOUNDARY  FOR  VALUES  TO  BE  PLOTTED  IN  EACH  OF 

C  THE  THREE  LEVELS  (AND  THE  PSEUDO  FOURTH  LEVEL) 

C  BOUNDU  -  The  upper  boundary  for  values  TO  BE  PLOTTED  IN  EACH  OF 

C  THE  THREE  LEVELS   (AND  THE  PSEUDO  FOURTH  LEVELl 

C  LETTER  -  THE  ONE  DIGIT  SYMBOL  TO  BE  PLOTTED  FOR  EACH  VARIABLE. 

C  J»LL  VARIABLES  FOR  ALTERNATIVES  ARE  PLOTTED  AS  -A-  AND 

C  ARE  PLOTTED  IN  THE  SAME  LEVEL  AS  THEIR  NORMAL 

C  COUNTERPART  TO  IDENTIFY  TflEM 

C  POINT  -  THE  ARRAY  WHICH  REPRESENTS  ONE  LINE  ON  THE  PLOT.      IT  IS 

C  DIVIDED  INTO  THREE   LEVELS   (INDEPENDENT  PLOTS),   WITH  A 

C  BASE   LINE  PRINTING  FOR  EACH  AT   POSITIONS   I,    42.  83 

C  AND  124.     THIS  LEAVES  40  POSITIONS  BETWEEN  THE  LINES 

C  FOR  PLOTTING  PURPOSES 

C  RAISE  -  THE  QUANTITY  NEEDED  TO  RAISE  THE  CURVE  TO  THE  PROPER  LEVEL 

C  SCALE  -   THE    SCALING  FACTOR   FOR  EACH  OF  THE  VARIABLES.  EACH 

C  INCLUDES  THE  NORMAL   SCALING  FACTOR.   20.0.  AND  ANY 

C  OTHER  FACTOR  DEEMED  NECESSARY.   AS  EXPLAINED  BELOW 

C  TOPLOT  -  An  array  of  level   INDICATORS  FOR  EACH  VARIABLE  TO  BE 

C  PLOTTED.      IF  IT   IS  ZERO.  THE  VARIABLE  WILL  NOT  BE 

C  PLOTTED.     THE  VALUE  OF  TOPLOT  MUST  BE   1,  2  OR  3  FOR 

C  ALL  VARIABLES  EXCEPT  STORAGE.     SINCE   STORAGE   IS  A 

C  NEGATIVE  VALUE  AND  PRINTS  BELOW  THE  BASE  LINE,    IT  MAY 

C  NOT  BE  PRINTED  AS  PART  OF  LEVEL   1.     IT  HAY,  HOWEVER. 

C  BE  ASSIGNED  TO  THE  PSEUDO  LEVEL  4.   AND  THUS  WILL 

C  PRINT  BENEATH  THE   TOP  HOST  BASE  LINE 

C 

COMMON  AIRTEMC (25.6) 

COMMON  CALDEF(25) .CDMAX(25) ,COVDEN( 251 
COMMON  DREADY(25I 

COMMON  ENGBAL ( 25 )  ,  E T ,E TDA I LY ( 25 . 1 2  I 

COMMON  FREEWAT(25) 

COMMON  IS0THRH(25,20) 

COMMON  LASTUSD(25) .LEVELl, LEVEL2 

COMMON  MMOD 

COMMON  NDAYSN0(25) , NOI V SBL , NSUB. NYE ARS 
COMMON  DNTREES(25) 

COMMON  PEAKPPT I  25,201 .PEAKWE ( 25 ,20 ) , PHASE ( 25  I , POTENT ( 24  I . 
1  PREWEQV(25) 
COMMON  RECHRG(25) 

COMMON  SIHTEMl (25.3) .SUBID(25,6) .SLPASP(25,24) 

COMMON  TCOEFF ( 25) , THRSHLO ( 25 ) , TOPLOT ( I 1 ) 

COMMON  VEGTYPE(25) 

COMMON  WEIGHTI25),WSHEDIDI6I 

COMMON  YEARS(20) ,YYMMDO 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDID 

INTEGER  YEARS, YYMMDD 

CQMMON/FORDATA/  F00TN0T(26) ,MAXLINE 
INTEGER  FOOTNOT 

COMMDN/M AS TER /DATE(3) ,THXMSTR,THNHSTR, PPTMSTR, PPTONOW, OBSHYDR, 
1   POTRAD.MSTREOF. lYR 
INTEGER  DATE 
COMMON/PLOT S/PLOT( 11) 

COMMON/ WATRBAL/ETFROM. EVAPOTR, GENRO,PREC I P,RADIN,RADLWN,RADSWN, 
1   TEMPMAX, TEMPMIN. WATERIN 
INTEGER   BOUNDL (4) ,B0UN0U(4) .POINT (124  I 
DIMENSION  LETTER(ll),RAISE(4),SCALE(ll) 
EQUIVALENCE   ( DATE! 2) ,IDAY) 

DATA  BOUNDL, BOUNDU, RAISE /I ,42,83,83,42,83,124, 124, 1.5,42.5,83.5, 
I  124.5/ 

DATA  P0INT/IH.,40*1H   ,1H.,40«IH   ,1H.,40*IH  ,1H./ 

C  HYDROGRAPH  -  MAX  VALUE  =  0.5,   MULTIPLY  BY  4  AS  WELL  AS  20 

DATA   LETTERI 1 ) .SCALE (1 ) /IHH.80.0/ 
C  WATER  EQUIVALENT  -  MAX  VALUE  =  30.0,   DIVIDE  BY   15  (1.33  -  20/15) 

DATA  LETTERI2),LETTER(7),SCALE(2),SCALE(7)/1HW,1HA,2»1.33/ 
C  INPUT  -   MAX   VALUE   =   2.0.   NO  EXTRA   SCALING  NEEDED 

DATA  LETTER(3I .LETTER(8) .SCALE (3) . SC AL E ( 8 ) / IHI . IHA. 2*20 . 0/ 
Q  EVAPOTRANSPIRATION  -  MAX  VALUE   =  0.5.    MULTIPLY   BY  4  AND  20 

DATA   LETTER(4),LETTERI9).SCALE(4),SCALE(9)/1HE, IHA, 2» 80.0/ 

C  STORAGE  REQUIREMENTS  -  MIN  VALUE  =  -5.3   (7.547  =  20/(5.3/2)1 

DATA   LETTER(5).LETTER(10).SCALE(5),SCALE(10)/IHS, IHA, 2*7. 547/ 

C  RUNOFF  -  MAX  VALUE  =  2.0,   NO  EXTRA  SCALING  NEEDED 

DATA  LETTER(6) .LETTER! 11  ) . SC ALE ( 6  I . SC ALE ( 1 1 ) / IHR . 1  HA . 2*2 0. 0/ 

C  SCALE  EACH  VARIABLE   THAT  IS  TO  BE  PLOTTED,   RAISE  IT  TO  THE  PROPER 

C   LEVEL  AND  IF   IT   IS  WITHIN  THE  BOUNDARIES  FOR  THAT  LEVEL.  STORE 

C   THE  CHARACTER   FOR  THE  PLOT 

00  20   I   =  l.U 
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IF( TOPLOT( I ) (  20,20.10 
10  J  =  TOPLOTII) 

IPOINT  =   (PLOT(I)   t  SCALE(I))   ♦  R*ISE(J) 
IFI [ .E0.5.0R. 1 .EQ. 10)   J   =   J  -  1 

IF( IPOINT. GT.BOUNDL( JI.AND.IPOINT.lt. 80UNDU1JI  I    PO  I  NT (  I PO  INT  1  = 
1  LETTER(I) 
20  CONTINUE 

C  WRITE  THE  DATE  BY  TENS 

IF( IDAV.E0.lO.OR.IDAY.EQ.2O.OR.IDAY.EQ.30l   GO  TO  30 
WRITE    111,910)  POINT 
910  F0RMATi9X124Al) 
GO  TO  40 
30   WRITE    111,9201   DATE, POINT 
920  F0RMAT(1X3I2,2H..12*A1) 
"lO  DO  50   I   =  2,123 
50  POINT! I  I   »  IH 
P01NT1<,2I   =  IH. 
P0INT(83)   =  IH. 
RETURN 
END 


Subroutine  RADCOMP 


SUBROUTINE  RADCOMP 

C  COMPUTE  THE   POTENTIAL  RADIATION   AT   THE   BASE  STATION 

C  

C  DICTIONARY 

C 

C  BETWEEN  -   THE   TOTAL   NUMBER  OF   DAYS   BETWEEN   THE  RESPECTIVE 

C  LOCATIONS  OF   THE   DATES   IN  -NDATE- 

C  DAYS  -   THE   NUMBER   OF   DAYS   THAT   HAVE   PASSED  SINCE   THE   BASE  DATE 

C  l-NOATE( ISUB)-)   TO   BE   USED   IN   THE  INTERPOLATION 

C  FRACTON  -  THE   FRACTIONAL   PART   NEEDED   IN   THE    INTERPOLATION  BETWEEN 

C  TABLE   VALUES   IN  THE  COMPUTATION  OF   THE  RADIATION 

C  ISUB  -   THE   SUBSCRIPT   OF   THE   BASE  TABLE  VALUE   USED  TO  OBTAIN  THE 

C  RADIATION   FROM  THE   TABLES  BY  INTERPOLATION 

C  NDATE  -  THE  DATES  OF  THE  TABLES  USED  IN  COMPUTING  THE  RADIATION 

C 

COMMON   AIRTEMC 125,61 

COMMON  CAL0EF(25 ) .CDMAXI 25) ,COVDENl 25  I 
COMMON  DREADYI25) 

COMMON  EnGBAl(25) ,ET,ETDAILYI25,12) 

COMMON  FREEMAT(25I 

COMMON  1SOTHRM(25.20) 

COMMON  LASTUS0(25) ,LEvEL1,LEVEL2 

COMMON  MHDD 

COMMON   NDAYSN0I25 ) .NDIVSBL.NSUB.NYEARS 
COMMON  0NTREESI25) 

COMMON  PEAKPPTI25,20)  ,  PE  AKWE  ( 25  , 20 1  .PHASE  I  25  ),  POTENT  (  2'.) , 
I  PREWE0VI25) 
COMMON  RECHRG(25) 

COMMON  SIMTEM1(25.3) . SUB  I D ( 25 , 6 ) . SL PA S P I  25 , 24 ) 

COMMON  TCOEFF(25l , THRSHLD 1 25 ) .TOPLOT ( 11 ) 

COMMON  VEGTYPEI25I 

COMMON  WEIGHTI25),WSHEDIDI6) 

COMMON  YEARS(20) ,YYMMOO 

INTEGER  DREAOY 

INTEGER  PHASE 

INTEGER  SUBIO 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  HSHEDIO 

INTEGER  YEARS, YYMMDD 

COMMON /MASrEK /DATE  131 , TMXMSTR , 1 MNMS TR , PPTMS T R . PPTONOW , OBSHY OR , 
I   POTRAD,HSTREOf , I YR 
INTEGER  DATE 
COMMON/RAD/FRACTON, I  SUB 

COHM0N/WATRBAL/ETFROM,EVAP0TR,GENRO.PRECIP,RADIN,RA0LHN,RADSHN, 
1  TEMPHAX.TEMPMIN.WATERIN 
DIMENSION  BETWEEN  I  2<.).  NDATE  1  24) 

DATA  BETWEEN/13..   15..   13..   15.,   14..    14..    15..    14.,    15.,  14., 

1  21.,    20..    15..    14..    15..    15..    14..    15..    14.,    14..    14..    14..  19., 

2  19./ 

DATA  NDATE/  110.  1 2 3 . 207 . 22 0 . 307 , 32 1 . 404 . 41 9 . 503 , 5 1  6 ,60 1  , 52 2 , 7 12 , 
I  727, 810, 625, 909, 923.1006,  1022. 11 05, II 19.  1 203.  1222/ 

C  PLACE   THIS   DATE   WITH  RESPECT   TO  THE  TABLES 

DO  10   I  =  1,24 
IFCNDATEd)  -  MMDDI  10,80,20 
10  CONTINUE 

t  4   NORMAL  TERMINATION  OF   THE  DO  LOOP  MEANS   THAT   THIS   DATE  FALLS 

C   BETWEEN  12/23  AND  12/31,    INCLUSIVE.     USING  THE  ARRAY   IN  CIRCULAR 

C   FASHION.    I/IO  (SUBSCRIPT   1)    IS  THE  CONTROLLING  DATE 

I   =  I 

GO  TO  30 

C  THIS  DATE  FALLS  BETWEEN  THE  ONES  AT  LOCATIONS   I   AND  I-l.     IF   1  IS 

C   1.  USE  24  FOR   I-l  SINCE  THE  ARRAY  IS  CIRCULAR 

20  ISUB  =  I  -  1 

IFIISUBI  30.30.40 

30  ISUB  =  24 

C  OBTAIN   THE    INTERPOLATION  FRACTION.      START   BY   DETERMINING  IF 

C   THIS  DATE   FALLS   IN  THE   SAME   MONTH   AS   THAT   AT   LOCATION    I   OR  I-l 

40  IFIDATEIl)   -   (NDATEUSUB)/100)  )  60.50.60 

C  IT   IS  THE  SAME  AS   I-l  AND  IT   IS  LARGER.   SO  SUBTRACT  THE   I-l  DATE 

C   TO  OBTAIN   THE   NUMBER  OF   DAYS   To   BE  USED  FOR  INTERPOLATING 

50  DAYS  -   MMDD  -  NDATEdSUBI 
GO  TO  70 

Q  ,T   IS  THE   SAME  AS   I,   BUT   IT   IS  SMALLER,   SO  SUBTRACT   IT  FROM  THE  I 

C   DATE.      THEN   SUBTRACT   THE   RESULT   FROM   THE   DAYS   BETWEEN    I   AND  I-l 

C   TO  OBTAIN   THE   NUMBER  OF   DAYS  TO  BE   USED  FOR  INTERPOLATING 

60   DAYS    =    NDATE(I)    -  MMDD 

DAYS  =  BETWEENdSUBI  -  DAYS 

C  COMPUTE   THE   INTERPOLATION  FRACTION 

70  FRACTON  =  DA YS/BE T WEEN  I  I  SUB  I 

POTRAD  -   POTENT! ISUB)   ♦   (POTENT! I)   -  POTENT ( I  SUB ) )   *  FRACTON 


RETURN 

C  THIS  DATE   IS   IN  THE  TABLE  -  NO  INTERPOLATION   IS  NECESSARY 

80  FRACTON   =  0.0 
ISUB   =  I 

POTRAD   =  POTENT!  I  ) 

RETURN 

END 

Subroutine  RDMSTR 

SUBROUTINE  RDMSTR 

C  READ   A  CARD  FROM  THE   MASTER  DECK 

COMMON  AIRTEMC(25,6I 

COMMON   CALOEF ( 25) ,CDMAX!25) ,COVDEN! 25) 
COMMON  DREADY(25I 

COMMON  ENGBAL125I , E T , ET DA  I LY ! 25 , 1 2 ) 

COMMON  FREEWAT(25] 

COMMON   I SOTHRMI25,20) 

COMMON  L4STUS0(25) ,LEVEL1.LEVEL2 

COMMON  MMDD 

COMMON  NDAYSNO! 25) , NO  I V SBL , NSU8 . NYE ARS 
COMMON  ONTREES(25) 

COMMON  PEAKPPT(25,2O),PEAKWE(25,20),PHASEI25),P0TENT(24l. 
1  PREWEQV(25> 
COMMON  RECHRG(25) 

COMMON   SIHTEM1(25.3).SUBID(25.6) . SLPASP ( 25 , 24 ) 

COMMON  TCOEFF (25) , THRSHLD ( 25 ), TOPLOT ( 11 ) 

COMMON  VEGTYPE(25) 

COMMON   WE IGHTI 25) .HSHEDIO (61 

COMMON    YEARS(20) , YYMMDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBin 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDIO 

INTEGER    YEARS, YYMMDD 

COMMON /MASTER /DATE (3) . TMXMSTR , TMNMS TR , PPTMST R , PPTONOW . OBSHYOR . 
1   POTRAD, MSTR6UF , I YR 
INTEGER  DATE 

C0MMON/WArRBAL/ETFROM,EVAP0TR,GENRO,PRECIP,RADIN,RADLWN,RADSWN, 
1  TEMPMAX,TEMPHIN,HATERIN 
DIMENSION  INI8) 
EQUIVALENCE    ( DATE  I  I ) , I N ( 1 ) ) 
READ    (5,9101  IN 

C  THE  FORMAT  DECLARATORS  ARE  TYPED  FOR  THE  EQUIVALENCED  WORDS  OF 

Q   RATHER  THAN  FOR  THE   INTEGER  ARRAY  ITSELF 

910   FORMAT I3I2,11X2F4.1.19X3F5.2) 
IF!E0FI5)I  20,20,10 
10  MSTREOF   =  1 
RETURN 

20   MMDD  =    (DATEIl)   »    100)    ♦  DATE(2I 

YYMMDD  =   (DATE(3)   *  lOOOOl   *  MMDD 

C  REDUCE   -PPTONOW-   TO   ITS  VALUE   BEFORE   THE  PRECIP   FOR   THIS  DAY  WAS 

C  ADDED  IN 

PPTONOW   =  PPTONOW  -  PPTMSTR 

C  SINCE   RADIATION  MEASUREMENTS   ARE  NOT   AVAILABLE,    COMPUTE   THE  BASE 

C   STATION  VALUE  HERE 

CALL  RADCOMP 

RETURN 

END 


Program  SELECT 


OVERLAY  lOLAYS.l.O] 
PROGRAM  SELECT 

C  SELECT   THE   METHOD   OF  OUTPUT   AND  READ  THE  WATERSHED  PARAMETERS 

COMMON   AIR TEMC ( 25. 5 ) 

COMMON  CALDEF(25) ,CDMAX(25) .C0VDENI25) 
COMMON  DREA0Y(25) 

COMMON   ENGBAL( 25) .ET,ETDAILY(25,12) 

COMMON  FREEWAT(25I 

COMMON   I SOTHRH(25,20I 

COMMON   LASTUSD(25) ,LEVEL1.LEVEL2 

COMMON  MMDD 

COMMON  NDAYSNO (25) , ND I  V SBL . NSUB . NYE ARS 
COMMON  0NTREES(25I 

COMMON   PEAKPPT(25.20) , PE AKWE ( 25 , 20 ) , PHASE  I  2 5), POTENT  I  24), 
1  PREWEQV(25) 
COMMON  RECHRG(25) 

COMMON  SIMTEMII25,3),SUBIDI25,6I.SLPASP(25,24) 

COMMON   TCOEFF I  25)  ,THRSHLD( 25 ), TOPLOT I  11 ) 

COMMON  VEGTYPE(25) 

COMMON  WEIGHT(25),WSHEDID!6) 

COMMON    YEARS(201 , YYMMDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDIO 

INTEGER   YEARS, YYMMDD 

COMMON /MASTER /DATE  131 , TMXMSTR , T MNMS TR , PPTMSTR , PP TONOW , OBSHYDR , 
I   POTRAD, MSTREOF, lYR 
INTEGER  DATE 

C 0MMON/WATRBAL/ETFROM,EVAPOTR,GENRO. PRECIP, RAD  I N.RADLWN.RAOSWN. 
1  TEMPMAX.TEMPMIN.WATERIN 
DATA  IERR/1/ 

C  READ  THE  OUTPUT  SELECTION  AND  CHECK  FOR  ERRORS 

READ   15.910)    I .LI ,L2.NDIVSBL 
910   F0RMAT(A6,1XI1 .1X11.1X12) 
IF ( I .EQ.6HSELECT I   GO  TO  20 
10  WRITE    (6,9201  I.Ll,L2.NDIVSBL 
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120  FDRMATI'lTHE  SELECT  CARD   (FIRST  CARD   IN  DECK)   IS   INCORRECT.     EITHE  INTEGER  OREADY 

IR    IT   DOES  NOT   HAVE   THE   WORD  -SELECT-   IN  COLUMNS   1-6,   THE   OUTPUT   OP  INTEGER  PHASE 

2TI0N«/»    IS  NOT  2,    3  OR   1,   OR   THE  -DATE   DIVISIBLE   BY-  VALUE    IS   INVA  INTEGER  SUBID 

3LID«/5X*C0LUHNS   1-6  =  » A6/ 5 X»OUTPUT   OPTION   (COL   81   =* 1 2/ 5X» AL TERNA  INTEGER  TOPLOT 

<,TIVE   (COL  101   =«I2/5X»0ATE  DIVISIBLE  BY  (COL  11-121   =*I3/»-J0B  ABO  INTEGER  VEGTYPE 

5RTE0»1  INTEGER  WSHEOID 

CALL  EXIT  INTEGER  YEARS. YYHMDD 

20  IF(L1.EQ.3.0R.L1.EQ.'.)   GO  TO  50  COHHON/CMPS I TO/COMPS  (  16 ) ,  YftTOT  ( 5  I 

IF(L1  -  21   10,30,10  COMHON/HASTER/OATEOl  ,THXHSTR,THNMSTR,PPTMSTR,PPTONOW,OBSHYDR, 

C  OPTION  2  -   CHECK   THE   DIVISIBILITY  1   POTRAD , MSTRE OF , I YR 

30  IF(NDIVSBL1   IClCiO  INTEGER  DATE 

AO   IFINDIVSBL  -   311    50,50,10  COMMON /WATR BAL /E TFROH , E VAPOTR , GENRO , PRECI P , RAO  I N , R AOL WN,  RAOSWN, 

50  LEVELl   =   LI  '    TEMPHAX , TEMPM I N , WATER  I N 

LEVEL2   =  L2  ♦   1  INTEGER  PEAKOAT 

C  READ  THE  PLOT  CONTROL  CARD  AND  CHECK  FOR  ERRORS  DATA  IERR,N/2«1/ 

READ   (5,9301   I, TOPLOT  C  READ  THE  WATERSHED  TITLE  AND  NUMBER  OF  SUBSTATIONS 

930  F0RMAT(A5,3Xll( IXIU  I  READ   (5,9001   WSMEO 1 0 ,NSUB , NYEARS 

IF( I .EQ.5HPLaTSl   GO   TO   60  900   FOR  HAT ( 6A 1 0 , 1 *X2 1  3 1 

WRITE   (6.9401  IF(NSUB.GE.1.AND.NSUB.LE.25I  GO  TO  10 

940  F0RKAT(»1THE   SECOND  CARD   IN  THE   DATA  DECK    IS   HOT   THE  PLOTS  CARD  -  WRITE    (6,9011  lERR.NSUB 

IJOB  ABORTEOtl  901   FORMAT ( I  1 , I  3 ,»   SUBSTATIONS   MAY  NOT   BE  RUN.     COLUMNS   76-77  OF   THE  W 

CALL  EXIT  lATERSHEO  10  CARD  MUST  BE  1  TO  25,   INCLUSIVE   (OR  THE  PROGRAM  MAY  BE 

C  HYDROGRAPH  2  REVISEDl*! 

60   IF(T0PL0T(11   -   31    BO, 80, 70  lERR   =  0 

70  I   =   lOHHYDROGRAPH  10  I F ( NYE ARS . GE . I . ANO. NYEARS . LE . 20  I  GO  TO  20 

WRITE    (6,9501    lERR, I ,T0PL0T( 1 1  WRITE    (6,9021    lERR, NYEARS 

950  FORMAT ( I  1 ,«THE   "AlO,*  MAY  NOT  BE  PRINTED  ON  BASE  LINE*I2,»  OF  THE  goj  FORMAT ( 1 1 , 1 3 , *  YEARS  MAY  NOT  BE  RUN.     COLUMNS  79-80  OF  THE  WATERSH 

1   PLOTS'l  lED  ID  CARD  MUST  8E  1  TO  20,   INCLUSIVE   (OR  THE  PROGRAM  MAY  BE  REVIS 

lERR  =  0  2EDI*I 

C  WATER   EQUIVALENT  ,ERR   =  o 

80   IF(T0PL0T(2I   -   31    100,100,90  C  READ   THE   POTENTIAL  RADIATION  CARDS    (2  CARDSl 

90   I    =   lOHWATER   EOV  20  READ   (5,9101   NAME. POTENT 

WRITE   16,9501    lERR, I .TOPLOTI 2  I  910  FORMAT ( A 10 , lOX 12F5. 0/20X12F5.0 1 

lERR   =   0  IFINAME.EO.IOHPOTENTIAL    1   GO  TO  30 

C  INPUT  WRITE    (6.9111    I  ERR 

100  IF(T0PL0T(31  -  31   120.120,110  911  F ORMAT ( 1 1 , » THE  POTENTIAL  RADIATION  CARDS  D0> NOT  FOLLOW  THE  WATERSH 

110   I   =    lOHINPUT  lEO   10  CARD«1 

WRITE    (6,9501    lERR.I .T0PL0T(31  lERR   =  0 

lERR  =  0  C  SEE   IF   THERE  ARE  SUBSTATIONS  TO  BE  READ 

Q  EVAPOTRANSPIRATION  30   IFINSUBl  210.210.60 

120   IFITOPLOTCI   -   31    1<.0.110.  130  60  NSETS   =  NSUB 

no  I   =  lOHEVAPOTR  C  READ  A   SET  OF   SUBSTATION  DESCRIPTORS.     START  WITH  THE   ID  AND 

WRITE    (6.9501    lERR. I .TOPLOT CI  C   CONSTANT  PARAMETERS 

lERR   =   0  70   READ   (5.9301    ( SUB  I  0 ( N, 1 ) . I  =  1 . 6 1 . TCOEFF ( N ) , COVDEN ( N I .CDMAX ( N I , 

C  STORAGE   REQUIREMENTS   (SINCE   THIS   IS   A  NEGATIVE   QUANTITY.    IT   MAY  i   WE  I GHT ( N I , THRSHLD ( N 1 , VEGTYPE ( N 1 

C  HOT   BE   PRINTED   AS   PART   OF   LEVEL   1,    BUT   MAY   BE   PART  OF  A  PSEUDO  930  FORMAT  I  6A 1 0 ,  IF 'i .  2  ,  F  2 .  0 .  1  X I  1 1 

t  LEVEL,   LEVEL  41  I F ( CDMA X ( N 1 . GE . C OVOE N( N ) 1   GO  TO  80 

140   IF(T0PL0T(51.LT.<..AND.T0PL0T(51.NE.ll   GO  TO  150  WRITE   (6,931)    I  ERR,  (  SUBI 0  (  N. I  1 . 1  =  1, 6 1 , COVDEN I N 1 .CDMAX ( Nl 

I   =   lOHSTORAGE  931  FORMAT ( 1 1 . «0N  THE   SUBSTATION  ID  CARD  ENTITLED  •6AlO/»  THE  COVER  DE 

WRITE   (6,9501    IERR,I ,T0PL0T(51  INSITY  SPECIFIED  IN  COLUMNS  65-58   (*F5.2,«I    IS  GREATER  THAN  THE  MAX 

lERR  =  0  2IMUM  COVER  DENSITY   IN  COLUMNS  69-72  (•F5.2,*l»l 

C  RUNOFF  lERR  =  0 

150  IF(T0PL0T(61  -  31    180,180,160  80  I F ( WE  I GHT( N 1 .GT.O. 0. ANO. WE IGHT ( Nl . LE . 1 . 0)   GO  TO  90 

160   I    =   lOHRUNOFF  WRITE    (6.9321    I E RR . WE  I GHT I N 1  . ( SUB  1 0 ( N, I  1 . 1= 1 . 5  I 

WRITE   (5,950)    lERR, I ,T0PL0T(61  932  F ORMAT ( 1 1 , • INVAL 1 0  WEIGHT   (•F5.2,»l    IN  COL  73-76  OF  SUBSTATION  10 

170  CALL  EXIT  ICARD  ENTITLED  *6AlO/»  WEIGHT  MUST  BE  BETWEEN  0.001  AND  1.0.  INCLUS 

180   IF(IERR)    170.170.190  2IVE*1 

C  BE    SURE   THAT  ALTERNATIVES   PRINT   ON  THE   SAME  LEVEL   AS  THEIR  lERR   =  0 

C             COUNTERPARTS  90  I F ( VEGTYPE ( N) .EO. 1 . OR. VEGTYPE ( N 1 .EQ. 2 1   GO  TO  100 

190  00  210   I   =   7.11  WRITE   (6.9331    I  ERR .  VEGTYPE  (  N  1 .  (  SUB  1 0  ( N .  I  ) .  I   1 . 6  1 

IF( TOPLOT! I  1 1   210.210.200  933  FORMAT ( I  1 , » I NVAL 10  VEG  TYPE   (*I1.*1    IN  COLUMN  80  OF  SUBSTATION  ID 

200  TOPLOTdl   =  T0PL0T(I-51  ICARD  ENTITLED  »5A10/*  VEGETATION  TYPE  =  1   (LODGEPOLE  PINEl,   =  2  (S 

210  CONTINUE  2PRUCE  FIR)*) 

C  READ  THE  WATERSHED  DESCRIPTORS  AND  PARAMETERS  IErr  =  0 

CALL  ROPARAM  C  READ   THE   INITIAL  CONDITIONS  CARD 

C  WRITE  THE  FIRST  LINES  OF  THE  PLOT  100  READ   (5,9401   NAME . S 1 MTE Ml ( N. 2 ) . PREW EQV ( N I . RECHRGI N ) . S IMTEM I ( N . 1 ) , 

REWIND   11  1  DREAOY(N) 

WRITE    (11.9601    WSHEOID. NSUB  940  FORMAT ( A  10, 10X4F5. 2 . I  5 1 

960  F0RHAT(»1»54X*WATER  BALANCE   SIMULATION*/  I F ( N AME . EO . lOH I Nl T I AL  COI   GO  TO  110 

1  1X6A10.41X*C0HP0SITE  0F*I3.*   SUBSTATIONS*//  WRITE    (6.9411    1 E RR , ( SUB  I  0 ( N , I ) . 1= 1 . 6 1 

2  •0»53X*LEGEND*/  941   FORMATdl.*  THE   INITIAL  CONDITIONS  CARD  DOES  NOT  FOLLOW  THE  SUBSTA 

3  1X2(12X«CHAR     0EFINITI0N«15X»RANGE«12X1/  ITION   ID  CARD  EN T I TL E 0* / 1 X6 A  1 0 ) 

4  14X«A         ALTERNATIVE   RESULT   FOR   NEAREST   CHARACTER  •                                      lERR   =  0 

5  lOXtR         RUNOFF-SIMULATED  GEN     0  TO  2.0   INCHES  •/                                   GO   TO  130 

6  14X.E         EVAPOTRANSPIRATION         0  TO  .2   INCH  *  C  CONVERT  THE  PACK  TEMPERATURE  TO  CALORIE  DEFICIT   (AS  A  POSITIVE 

7  10X«S         STORAGE  REQUIREMENTS     -5.3  TO  0   INCHES  */  C  QUANTITY),   INCLUDE   THE  WEIGHTED  RECHARGE  REQUIREMENT  FOR  USE  IN 

8  14X«H         HYDROGRAPH-OBS  GENRO     0   TO   .5  AREA    INCHES         •  C   THE   -CHANGE    IN   STORAGE-  COMPUTATION.    AND  DEFINE  THE  GROUND 

9  10x»W         WATER  EQUIV  OF  PACK       0  TO  30.0  INCHES  */  C  TEMPERATURE  FOR  THE  SIMULATION  MODEL 

A   14X*I         INPUT  -  MELT  OR  RAIN     0  TO  2.0   INCHES  •/*0*)                 no  CALDEF(N1   =  -   SIHTEM1IN,21   «  PREWEOV(Nl   »  1.27 

C  RETURN   TO   THE   MAIN  OVERLAY   FOR  THE   LOADING  OF  THE  NEXT  OVERLAY  YRT0T(41   =   YRT0T(41    *    (RECHRG(Nl    *  WEIGHT(NI) 

END  SIMTEMIIN.3)   =  -  1.5 

C  READ   THE   DAILY  ET  VALUES 

O     1                         l-»T-»T->  A  Tl  *  n/T  130  READ   (  5.950)   N  AME  ,  (  E  TOA I  L  Y  (  N .  I  I  .  I  =  1  ,  1 2 ) 

bubroutme  RDPARAM  950  format(aio,ioxi2f5.4i 

IFINAME.EQ. lOHDAI LY  ET      )   GO  TO  140 

SUBROUTINE   ROPARAM  VRne    (5,951)    1  ER  R  ,  (  SU  B  1 0  (  N  ,  1)  ,  I  =  1 ,  6  I 

C  READ   THE   WATERSHED   PARAMETERS   AND  DESCRIPTORS  951   FORMATdl,*   THE   DAILY  ET  VALUES  CARD  DOES  NOT  FOLLOW  THE   INITIAL  C 

COMMON  AIRTEMC(25,61  lONOITIONS  CARD   IN  THE  CARDS  FOLLOWING  THE  SUBSTATION   ID  CARD  ENTIT 

COMMON  CALOEFI 25) ,C0HAX(25) , COVDEN!  25)  21 E 0* / I X6 A  1 0 ) 

COMMON  0REA0Y(25)  ,51,1,  ^  q 

COMMON  ENGBAL(251,ET,ETDAILY(25,121   rjad  ThE  AIR  TEMPERATURE  COEFFICIENTS 

COMMON  FREEWAT(251  l^n  READ    (5,960)   NAME , ( A  I RT EMC ( N, I ) . I  =  1 . 5 1 

COMMON   IS0THRM(25.20)  ,(,0  FORM  AT  |  A  1 0  .  10X5F  5.  3  ) 

COMMON  LASTUS0(251. LEVELl. LEVEL2  I F ( NAME . EQ . lOHA I R   TEMP  CI   GO  TO  150 

COMMON  MMOD  WRITE    (5.9611    I  ERR . ( SUB  1 0 1 N , I  I , I  =  I , 5 ) 

COMMON  NOAYSNO(25).N01VSBL. NSUB. NYEARS  ,61   FORM A T ( I  1 , « THE   AIR   TEMPERATURE  COEFFICIENTS   CARD  DOES  NOT   FOLLOW  T 

COMMON  ONTREES(25)  IHE  DAILY  ET  CARD  IN  THE  CARDS  FOLLOWING  THE  SUBSTATION  ID  CARD  ENT 

COMMON  PEAKPPT (25. 20 ) .PEAKWE (25.201  .PHASE ( 25 ). POTENT (  24 ) .  31 TLEO*/ 1 X5A 10 1 

1               PREWE0V(251  lERR  ,  q 

COMMON  RECHRG(25)  q  DEAD  THE  SLOPE/ASPECT  ADJUSTMENT  FACTORS 

COMMON  SIMTEM1(25.3).SUBID(25.6  1.SLPASP(25,241 

COMMON  TCOEFF (25) ,THRSHLD(25) , TOPLOT! 11  1  150  READ  (5,9701   NAME . ( SLPASP ( N. I ) , I  =  1 . 241 

COMMON  VEGTYPE(251  970  FORM AT ( A5 . 2 X24F 3 . 2 1 

COMMON  WEIGHT(251 ,WSHEDI0(6)  I  F ( NAME . EQ . 6HSL P/ AS  1   GO  TO  160 

COMMON  YEARS(201,YYMM00  WRITE    (6,9711    I  ERR . ( SUB  1 0 ( N. I  I . I  =  1 . 5 1 
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971  FORH&T( I1.»THE  SLOPE/ASPECT  CORRECTION  FACTORS  CARD  DOES  NOT  FOLLO 
IW  THE  AIR  TEMPERATURE  COEFFICIENTS  CARD  FOLLOWING  THE  SUBSTATION  I 
2D  CARD4/1X6A10I 
lERR  =  0 

C  READ  THE  PEAK  HATER  EQUIVALENT.    ITS  DATE  AND  THE  DATE  BY  WHICH  THE 

C  —   PACK  MUST   BE    ISOTHERMAL  FOR  EACH  YEAR 

160  00  190   I  =  l.NYEARS 

READ   (5,9801   PEAKWE ( N, I  I • TEMP , PE AKOAT . I  TEMP 

980  FORMAT! 20X2F5. 2. 1X16. 1X16) 

C  FOR   RAPID   AND   MORE   ACCURATE   BINARY   COMPARES    IN   SUBROUTINE  GENDATA. 

C   LOWER  THE   PEAK  PRECIP  SLIGHTLY 

PEAKPPTIN.n    =   TEMP  -  0.001 

IFIN  -   1 )   170, 170, 160 
170   YEARSII)   =  MOOIPEAKDAT.lOOl 

GO  TO  190 

180   IF( YEARS! I  1 .EQ.MOO(PEAKDAT, 100) 1   GO  TO  190 
K   '  HOD(PEAKDAT,100) 

WRITE    (6,9811    lERR , (SUBIDIN, J) . J=l,6) .YEARS  1  I ) .K. ( YEARS! J ) . J  =  l, 
1  NYEARS) 

981  FORMAT! I  I, •THE   PEAK   W.E.   CARDS  FOLLOWING   THE   SUBSTATION   ID  CARD  EN 
ITITLED  •6A10/»  DO  NOT  CORRESPOND  TO  THE  YEARS   AND/OR  ORDER  OF  THE 
2FIRST   SET  READ.      WHEN«I3,»  WAS   EXPE CT ED, • I  3 , *  WAS   READ.*/*   THE  FIR 
3ST  SET  IS0»20I'>) 

I  ERR  =  0 

C  CONVERT   THE   ISOTHERMAL  DATE   TO  YYHHDD   FORMAT   FOR   RAPID  COMPARES 

190   ISOTHRH(N,I)   -    (10000  *   ITEMP)    -   (999999  *  IITEHP/lOOl) 

C  GO  ON  TO  THE  NEXT  SET 

N  =  N  <•  lERR 
NSETS  =  NSETS  -  1 
IFINSETS)  200,200,70 
200  IF(IERR)  210,210,220 
210  WRITE  (6.990) 

990  FORHAT(»0«/»OJOB  ABORTED  FOR  ABOVE   REASONIS)   WHILE   OPERATING   IN  SU 
IBROUTINE  RDPARAM*) 

CALL  EXIT 

C  INITIALIZE    THE   VARIABLES   THAT  HAVE   NOT   BEEN  DEFINED  ABOVE 

220  DO  230  N  =  1.25 

ENGBALINI   =  -1.0 

FREEWAT(N)    =  0.0 

LASTUSDINI    =  0 

NDAYSNO(N)   =  0 

ONTREES(N)   =  0.0 
230   PHASE(N)   =  0 

C  SUM  THE  WEIGHTS  To  CERTIFY  THAT  A  COMPLETE  SET  OF  PARAMETERS  IS 

C   INCLUDED 

SUM  =  0.0 

DO  2'.0  N  =  l.NSUB 
240   SUM  =   SUM  *  WEIGHTIN) 

IF( SUM. GT.O. 97. AND.SUM.lt. 1.03)  RETURN 

WRITE    (6,991)  lERR.SUM 

991  F0RMAT(I1,*THE   TOTAL  OF   THE   WEIGHTING  FACTORS   («F6.3,«)    IS  OUTSIDE 
1   THE   TOLERABLE   TOTALS,   0.97   TO  1.03*) 

GO  TO  210 
END 


COMMON  FREEWATI25) 

COMMON  ISOTHHM(25,20) 

COMMON   LASTUSDI 25 ) ,LEVEL1,LEVEL2 

COMMON  MMDD 

COMMON  N0AYSN0I25) , NDI VS8L , NSUB , N YE ARS 
COMMON  0NTReES(25) 

COMMON   PEAKPPT(  25,  20  I  .PEAK  WE  I  25,  20  1, PHASE  (25), POTENT!  2'.  I, 
1  PREHE0VI251 
COMMON  RECHRGI25) 

COHMPN  SIMT6MU25,3),SU8ID(25,6),SLPASPI25,2'.) 

COMMON  TC0EFF!25), THRSHLOl 25) ,TOPLOT! 11 ) 

COMMON  VE&TYPE!25) 

COMMON    HEIGHT!25) ,WSHEDID!6) 

COMMON   YEARS  120 ) lYYMMDD 

INTEGER  OREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDID 

INTFGEB   YEARS, YYMMDD 

COMMON/CMPSITO/COMPS! 16) ,YRTOT(  5) 

COMMON/MASTER /DATE ( 3  I , TMXM S TR , TMNHS TR , PPTMS TR , PPTONOW , OBSHYOR, 
1    PaTRAO,MSTREOF . I YR 
INTEGER  DATE 
COMMON/ PLOTS/ PLOT (6 1 

COMMON /WATBBAL/BTFROM.EVAPOTR.GENRO, PRECIP, RAD  I N,RADLWN,RADSHN, 
1  TEHPMAX.TEMPMIN.WATERIN 

C  HAKE    THE   PASS   THROUGH   THE   WATER   BALANCE  ROUTINES 

CALL   WATBAL    ( C ALDEF I N 1 , COM AX ( N ) , COvOEN ! N ) , DREADY ( N ) , ENGBA L I N 1 . 

1  FREEWAT !N 1 ,LASTUSD(N) ,NDAYSNO(N) ,ONTREES(N) ,PHASE(N) ,PREWEOV(N) , 

2  RECHRGI N) ,SIMTEMl(N,l) , S I MT EM  I ! N , 2 ) , S I MT EMI  1 N, 3 ) , TCOEFF ! N ) , 

3  THRSHLD !N) , VEGTYPE I N) ) 

C  CHECK   THE   MANDATORY    ISOTHERMAL  DATE 

IF!YYMMDD  -    I  SO THRM I N , I YR )  1  20,10,20 

i;  TURN  OFF   THE   DIFFUSION   MODEL,    SET  THE   PACK   TO  0  C   AND  ADJUST  THE 

C   ENERGY  BALANCE  ACCORDINGLY 

10  DREADYIN)   =  -I 

FNCBALIN)   =  ENCBAL(N)   t  CALOEFIN) 

RADLWN   =  RADLWN  *  CALDEFIN) 

CALOEFIN)   =  0.0 

C  WEIGHT  AND  STORE  THIS  DATA 

20  WT   =  WEIGHTIN) 

COMPSll)    =  COMPSIl)    +    IPREHEQVINl    *  HT) 

C0MPS(2)   =  C0MPS(2)   ♦   IRECHRG!N)   •>  WT ) 

C0MPS!13)    =  C0MPSI13)    •■    IPRECIP   «  HT) 

COMPSd'')    =  COMPS(ll)    *    (WATERIN  «  WT ) 

COMPSIIS)  =  C0MPS(15I   ♦   (EVAPOTR  «  WT] 

C0MPS(161   =  C0MPS(161   *   (GENRO  «  WT) 

CALL   ETCODE    ( E TFROM , COMPS ( 6 ) ) 

RETURN 

END 


Subroutine  WRITE2 


Program  INTSUM 


OVERLAY  (OLAYS.2,0) 
PROGRAM  INTSUM 

C  INTERVAL   SUMS  OF   VALUES  CONCERNED  WITH  WATER  -   LOAD  APPROPRIATE 

C   OPERATING  PROGRAM  TO  WORK  WITH   THE   SUBROUTINES    IN   THIS  OVERLAY 

COMMON  AIRrEMC(25,6) 

COMMON  CAL0EF(25),CDMAX(25l,COV0EN( 25) 
COMMON  DREA0Y(25I 

COMMON  ENGBALI25) ,ET,ETDAILY(25,12) 

COMMON  FREEWAT(25) 

COMMON  IS0THRH(25,20) 

COMMON  LASTUSD(25), LEVEL), LEVEL2 

COMMON  MHOD 

COMMON  NDAYSN0(25) , NOi VSBL , NSUB , NYE ARS 
COMMON  ONTREES(25l 

COMMON  PEAKPPT(25,20) .PEAKWE (25,20) , PHASE ( 25 ), POTENT ( 2« ) , 
1  PREHEQV(25I 
COMMON  RECHRG(25) 

COMMON  SIMTEM1(25,31,SUBIDI25,6),SLPASP(25,2'.) 

COMMON  TCOEFF ( 25) , THRSHLD (25) ,TOPLDT( 11 ) 

COMMON  VEGTYPE (25) 

COMMON  HEIGHTI 25) ,HSHEDI0(6) 

COMMON  YEAR$(20), YYMMDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  HSHEDID 

INTEGER  YEARS, YYMMDD 

COMMON /MASTER /DATE (3) , TMXMSTR .TMNMSTR, PPT  MS TR, PPTONOW , OBSHYOR , 
1  POTRAO,MSTREOF,IYR 
INTEGER  DATE 

COMMON /WATRBAL/ETFROM, EVAPOTR, GENRO, PRECIP, RAD  IN, RADLWN, RADSWN, 
1  TEMPMAX,TEMPMIN,HATERIN 
CALL  OVERLAY   ( 5H0LAYS ,2 , LEVEL2 I 
C  RETURN  TO  THE   MAIN  OVERLAY   TO  TERMINATE  NORMALLY 

END 


Subroutine  NORM2 


SUBROUTINE  N0RM2  (Nl 

C  THE  VERSION  OF  THE  NORMAL  SIMULATION  SUMS  THOSE  VALUES  CONCERNED 

C   WITH  WATER  FOR  OUTPUT  ON  DATES  DIVISIBLE  BY  THE  SPECIFIED  NUMBER 

COMMON  AIRTEMCI25,6I 

COMMON  CALDEF (25) .COMAX ( 25 ) ,COVDEN( 25 ) 
COMMON  DREADY(25) 

COMMON  eNGBAL(25) ,ET,ETDA1LY(25, I?) 


SUBROUTINE  WRITE2   (K, COMPS) 

C  THIS   VERSION  OF   THE  OUTPUT   ROUTINE   PRINTS    INTERVAL   AND  YEAR-TO- 

C   DATE   SUMS   OF   ALL   VALUES  CONCERNED  WITH  WATER 

COMMON  4IRTEMC(25,6I 

COMMON   CALDEF (25), COMAX (25) ,COVDEN(25) 
COMMON  DREADY(25) 

COMMON   ENGBAL ( 25) ,ET.ETDAILY125, 12) 

COMMON  FREEWAT(251 

COMMON    I SDTHRH(25,20) 

COMMON  LASTUSD(25) ,LEVEL1,LEVEL2 

COMMON  MMOD 

COMMON  NDAYSNO( 25) , NDI V SBL , N SUB , NYE AR S 
COMMON  0NTREES(?5) 

COMMON  PEAKPPT  (25.20)  ,  P  E  AKWE  (  25  . 20  )  ,  PHASE!  25) ,  POTENT!  2^1 1, 
1  PREWEQV(25) 
COMMON  RECHRG(25) 

COMMON  SIM  IE  Ml (25, 3), SUB  ID (25, 6) .SLPAS P ( 25 , 24 ) 

COMMON   TCOEFF (25)  .THRSHLD! 25 ) , TOPLOT!  11) 

COMMON  VEGTYPE(25) 

COMMON   HEIGHT!25) .HSHE0IDI6) 

COMMON    YEARS  I  20 1  .YYMMDD 

INTEGER  OREADY 

INTEGER  PHASE 

INTEGER  SUBIO 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  HSHEOID 

INTEGER    YEARS, YYMMDD 

COMMON/FORDATA/   F00TN0T(26) .MAXLINE 
INTEGER  FOOTNOT 

COMMON /MASTER /DATE (3), TMXMSTR, TMNMSTR, PPTMSTR.PPTONOH, OBSHYOR, 
1   POTRAD,MSTREOF , lYR 
INTEGER  DATE 

COMMON / HATR8AL/ETFR0M, EVAPOTR, GENRO. PREC I P, RAD  IN, RADLHN, RADSWN, 
1  TEMPMAX.TEMPMIN.WATERIN 
DIMENSION  C0MPS(12) ,FR0M(7I 
DATA  LINES/-1/ 

DATA   FR0M/3HC      , 3H   S   ,3HCS   , 3H     E,3HC   E,3H  SE,3HCSE/ 
C0MPS(6)    =   FR0M!C0MPS(6) ) 

C  DETERMINE   HOW  TO  WRITE   THE  LINE 

IF !K.E8. IHt 1   GO  TO  10 
HRITE    !6,910)    K, COMPS 
910  F0RMAT!1XA10,2F10.2,6X2F8.2,FB.'.,1XA3,F9.2,7X2F8.2,F10.4, 1X2F10.2) 
LINES   =   LINES  -  1 
RETURN 

C  CHECK   THE   LINE  COUNTER 

10  IFILINESI  20,20.30 

20  HRITE    (6,920)   FOOTNOT , H SHED  I  0, N SUB 
920  F0RMAT!»0*I3A10/1X13A10/«1»5<.X»WATER  BALANCE  SIMULATION*/ 
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1  1X6410, <.1X»C0MP0SI  TE  0F*I3.«  SUBSTATIONS'/ 

2  «0«16xeC  U  R  R  E  N  T«12X»I   NTERVAL         TOTAL  S*13X 
3»---YEAR         TO         DATE-     -  -*/ 

A  lAXtSNOWPACK     RECHARGE»23X«EVAPOTRANS     GENE R ATED*3ex 

5  'GEN         CHANGE  IN«/ 

6  16X*H.   E.  REO  PRECIP        INPUT  FROM  RUNOFF 

7  PRECIP        INPUT     EVAPOTRANS         RUNOFF     RECHRG  R0»/) 
LINES   =  HAXLINE 

30  WRITE    (6,930)  DATE,COMPS 
930  FORMAT!*  *3I3,lX2Fl0.2,6X2Fe.2,F8.<.,lXA3,F9.2,7X2F8.2,F10.1,lX 
1  2F10.2) 
LINES  =  LINES  -  1 
RETURN 

C  DOUBLE  SPACE  BETWEEN  YEARS 

ENTRY  WRITOT 
WRITE  16,9101 
LINES   =   LINES  -  1 
RETURN 
END 


Program  INTSUMO 


OVERLAY  I0LAYS,2,1I 
PROGRAM  INTSUMO 

C  OPERATING  PROGRAM  FOR  PRINTING  INTERVAL  SUMS  FOR  THE  NORMAL 

C   SIMULATION   WITHOUT  ANY  ALTERNATIVES 

COMMON  AIRTEHC (25,61 

COMMON  CALDFF(251 ,CDMAX(25 ) ,COVOEN( 25) 
COMMON  0REA0Y(25) 

COMMON  ENGBAL(25  1 ,ET,ETDAILY(25,12) 

COMMON  FREEHAT(25) 

COMMON  IS0THRM(25,20) 

COMMON  LASTUSD(25) ,LEVEL1,LEVEL2 

COMMON  MMOD 

COMMON  NDAYSNO(251 , NDl VSBL , NSUB , NYEARS 
COMMON  0NTREES(25) 

COMMON   PEAKPPTI  25,201 ,PEAKWE (25,20) .PHASE  I  25) , POTENT (24). 
1  PREWEQV(25) 
COMMON  REChRG(25I 

COMMON   SIMTEM1(25,3),SUBI0(25,6  1,SLPASP(25,24) 

COMMON  TCOEFFI 251 ,THRSHLD(251 ,T0PL0T(11 ) 

COMMON  VEGTYPE(25) 

COMMON  WEIGHT(25) ,WSHE0ID(6) 

COMMON  YEARS(20l,YYMMD0 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTFGER  VEGTYPE 

INTEGER  WSHEDID 

INTEGER   YEARS, YYHHDD 

COMMON /CMP S I TO/COMPS ( 16 1 , YRTOT( 5  I 

COMMON/FORDATA/  F00TN0TI26) ,MAXLINE 

INTEGER  FOOTNOT 

COMMON /MASTER /DATE ( 3  I , T MXMSTR , T MNMSTR , P PTMSTR , PPTONOW ,08 SHYDR , 
1   POTRAD,MSTREOF, lYR 
INTEGER  DATE 
COMMON/PLOTS/ PLOT ( 11 ) 

COMMON/WATRBAL/ETFROM,EVAPOTR,GENRO,PRECIP,RAOlN,RADLWN,RADSWN, 
1  TEMPMAX,TEMPMIN,WATERIN 
INTEGER  F0aTNTE(26) 

DATA   F0OTNTE/26OHN0RMAL   SIMULATION  ONLY 

1 


C  DEFINE  THE  FOOTNOTE 

DO  1  N  =  1,26 
1  FOOTNOTINI   =  FOOTNTE(N) 
C  MAXLINE   =  52  -  NUMBER  OF  ALTERATIONS 

MAXLINE   =  52 
C  WRITE   THE   FOOTNOTE   ON  THE  PLOT 

WRITE   (11,900)  FOOTNOT 
900  FORMAT ( /*0*13A10/1X13A10//1 

C  READ   A   MASTER  CARD 

10  CALL  ROMSTR 

IFIMSTREOFl  20,20,170 

C  GENERATE   THE  DATA  AND  PERFORM  THE   SIMULATION  FOR  EACH  SUBSTATION 

20  DO   30  N  =  1,NSUB 

CALL  GENDATA  (N) 

CALL   N0RM2  (N) 
30  CONTINUE 
C  STORE   THE   INFORMATION  FOR  THE  PLOTS 

PLOT(l)   =  OBSHYDR 

PL0T(2)    =  COHPS( 1 ) 

PL0T(3)    -  COMPS(  l*,) 

PLOT!/,)   =  C0MPS(15) 

PL0T(5)   =  C0MPS(21 

PL0T(6)   =  COMPSI 16) 
C  ADO   THESE   VALUES   INTO   THE   INTERVAL  LOCATIONS 

C0MPS(3I   =  COMPSO)    ♦  C0HPS(13) 

C0MPS(4I   =   COMPSCil    +  COHPSIl*) 

CnHPS(5)   =  C0MPS(51   +  C0MPS(151 

C0MPS(7I    =   C0MPS(7)    +  C0HPS(16) 
C  ZERO  THE^AILY  COMPOSITE  LOCATIONS 

DO  40  N  =  13,16 
40  COHPS(N)   =  0.0 
C  SEE   IF   THIS  DAY   IS   TO  BE  PRINTED 

IFIMOD    (0ATE(2) ,N0IVSBL) )  70,50,70 

C  ADD  THESE   TOTALS   INTO  THE  YEAR-TO-OATE   LOCATIONS,   GET   THE  CHANGE 

Q   IN  jHE   RECHARGE  REQUIREMENTS   AND  PRINT   THE  LINE 

50  COMPSIS)    =  COMPSiai    *  C0MPS(3) 

C0MPS(9)   =  C0HPS(91   +  COMPSC) 

COMPSIIOI   =  COMPS(IO)   +  C0MPS(5) 


COHPSdll   =  COMPS(II)   t  C0MPS(7) 
C0HPS(12)    =  C0HPS(2)   -  YRTOTC) 
CALL  WRITE2  (1H$,C0MPS) 

C  ZERO  THE   INTERVAL  ACCUMULATING  LOCATIONS 

00  60  N   =  3,7 
60  COMPS(N)   =  0.0 

C  PERFORM  THE  PLOTS  BETWEEN  APRIL   1  AND  SEPTEMBER  30 

70  IFIMMDD  -  401)  160,140,60 
80   IF(930  -  MMDO)  160,90,150 

C  ON  9/30,   STORE   THE  CURRENT  RECHARGE  REQUIREMENT,   ZERO  THE  YEARLY 

C   ACCUMULATING  LOCATIONS,   AND  RESET  THE  DIFFUSION  MODEL  SWITCHES 

90  YRTOTIA)   =  C0MPS(2) 
00  100  N  =  8,11 
100  COHPSIN)   =  0.0 

DO   110  N  =  1,NSUB 
110  DREADY(N)   =  0 

C  CHECK   THE  YEARS  TO  BE  SURE  ALL  DECKS  STILL  CORRESPOND  TO  THE 

C  PARAMETER  DECK 

IF(DATE(31  -  YEARS(IYR))  120.130,120 
120  WRITE   (6,910)   DATE ( 3 ) . I YR . I YR , YEARS ( I YR ) 

910  F0RMAT(»1THE  DATA  DECKS   ARE   OUT   OF  PHASE  WITH  THE  WATERSHED  PARAHE 
ITER  DECK.     THE  WATER  YEAR  JUST  ENDING   (19*12.*)    IS  DECK  NUHBER»I3, 
2  *,»/«     BUT  SPECIFIED  CONDITIONS  CARD  NUMBER*I3,*  FOR  EACH  SUBSTAT 
3I0N   IS  FOR  19*12) 
CALL  EXIT 
130  lYR  =   lYR  ♦  I 

C  DOUBLE   SPACE  BETWEEN  YEARS 

CALL  WRITOT  (1H$.YRT0T) 
GO  TO  150 

C  ON  4/1,   WRITE  THE  ORDINATE  LINE 

140  WRITE   ( 11,920) 
920  F0RMAT(9X124(1H.  ) ) 
150  CALL  PLOTTER 

C  ZERO  THE  CURRENT  LOCATIONS 

160  COMPS( I)  =  0.0 
C0HPS(2)  =  0.0 
GO  TO  10 

C  ALL  CARDS  HAVE  BEEN  READ.      IF  A  LINE  WAS  NOT  PRINTED  ON  THE  LAST 

C   DAY  PROCESSED.  DO  IT  NOW 

170   IFIMOO   (MMDD.NDIVSBL I )  180.190.180 
180  COMPS(B)   =  C0MPS(8I   +  COMPSO) 
C0MPS(9)   =  C0HPS(91   +  C0MPS(4) 
COMPS(IO)   =  COMPS(IO)  +  C0MPS(5I 
COMPS(ll)  =  COMPS(ll)   +  C0MPS(7) 
C0MPS(12)   =  C0MPS(2)  -  YRT0T(4) 
CALL  WRITE2  (IHt.COMPS) 

C  WRITE  THE  FOOTNOTE  ON  THE  LAST  PAGE 

190   WRITE    (6,930)  FOOTNOT 
930  FORMAT(»0*13AlO/lX13A10) 

C  RETURN  TO  THE  MAIN  OVERLAY  FOR  NORMAL  TERMINATION 

END 


Program  DAILY 


OVERLAY  (0LAYS,3t0) 
PROGRAM  DAILY 

C  COMPOSITE  DAILY  OUTPUT  -  LOAD  APPROPRIATE  OPERATING  PROGRAM  TO 

C   WORK   WITH   THE   SUBROUTINES   IN   THIS  OVERLAY 

COMMON  AIRTEMC(25,6) 

COMMON  CAL0EF(251 ,CDMAX(25) .C0VDEN(25) 
COMMON  DREADY(251 

COMMON  EnGBALI 251 , ET . ETDA I LY ( 25 , 12 ) 

COMMON  FREEWAT(25) 

COMMON  IS0THRM(25,20l 

COMMON  LASTUSD(25) ,LEVEL1.LEVEL2 

COMMON  MMOD 

COMMON  N0AYSNO(25) ,NDI VSBL . NSUB. NYEARS 
COMMON  ONTREES(25) 

COMMON  PEAKPPT  (25,  20'),PEAKWE(25,20)  ,  PHASE  ( 25 ),  POTENT  (  24)  , 
1  PREWEQV(25) 
COMMON  RECHRG(25) 

COMMON  SIMTEM1(25,3),SUBID(25,6),SLPASP(25,24) 

COMMON  TCOEFF( 25) , THRSHLO( 25 ) ,TOPLOT( 11 ) 

COMMON  VEGTYPE(25) 

COMMON  WEIGHT(25) ,WSHE0IDI6) 

COMMON   YEARS(20) ,YYMMDO 

INTEGER  OREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEDID 

INTEGER   YEARS, YYHMDD 

COMMON/MAS TER/ DATE (3) >TMXMSTR,TMNMSTR,PPTMSTR.PPTONOW,OBSHYDR, 
1   POTRAD,MSTREOF, lYR 
INTEGER  DATE 

COMMON /W4TRBAL/ETFR0M.EVAP0TR.GENR0, PRECIP. RAD  IN. R AOL WN.RADSWN, 
1  TEMPHAX,TEMPMIN,WATERIN 
CALL   OVERLAY   ( 5H0LAYS , 3 ,LEVEL2 ) 
C  RETURN   TO  THE   MAIN  OVERLAY  To  TERMINATE  NORMALLY 

END 


Subroutine  NORMS 


SUBROUTINE   N0RM3  (N) 
-THIS   VERSION  OF   THE   NORMAL   SIMULATION  MAINTAINS   ALL  INFORMATION 

-  NECESSARY  OF   THE   PRINTING  OF  ONE  LINE   PER   DAY   FOR   THE  WATERSHED 

-  COMPOSITE 

COMMON  AIRTEMC(25,61 

COMMON  CALOEF (25) , CDMAX ( 25 ) .COVDEN 1 25 ) 

COMMON  DREADY(25) 

COMMON  ENGBAL (25 ) .ET,ETDAILY(25, 12) 
COMMON  FREEWAT(25) 
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COMMON  1SOTHRM(25,20I 

COMMON  L»STUSD(25l,LEVELl,LEUEL2 

COMMON  MMDD 

COMMON  NOAySNO(25l ,NDI VSBL ,NSUB ,NYE ARS 
COMMON  0NTREES(25) 

COMMON  PEAKPPT(  25,20),PE»KWE(25,20I  ,PHASEI25I  .P0TENTI2'(I  , 
1  PREHE0V(25) 
COMMON  RECHRS(25I 

COMMON  SIMTEMl 125,3) ,SUBID( 25,6) ,SLPASP(25t24l 

COMMON  TC0EFFI25) ,THRSHL0(25) ,TOPLOT( tl I 

COMMON  VEGTYPE125) 

COMMON  WEIGHT125) ,WSHE0I016) 

COMMON  YEARSI201,yyMHDO 

INTEGER  DREADV 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  HSHEOID 

INTEGER  YEARS, YYHMDO 

COMHGN/CMPSITO/COMPSl 16),YRT0T( 5) 

COMMON/MASTER /OATEI 3  I, TMXMSTR,TMNMSTR,PPTMSTR,PPTGNOW,OBSHYDR, 
1  POTRAD.MSTREOF,IYR 
INTEGER  DATE 
COMMON/ PLOTS/ PLOT  (6) 

COMHON/HATRBAL/ETFROH,EVAPOTR,GENRO,PRECIP,RA01N,RAOLHN,RADSHN, 
I  TEMPMAX,TEMPMIN,HATERIN 

C  MAKE   THE   PASS  THROUGH  THE  WATER   BALANCE  ROUTINES 

CALL   HATBAL   I C ALDE F 1 N ) , COM Ax ( N ) , COVDEN ( N ) , ORE ADY ( N 1 ,ENG8AL ( N ) , 

1  FREEWATIN ) ,LASTUSD(N) .NDAYSNOiN) .ONTREESIN) .PHASE  IN) ,PREWEQV(N I , 

2  RECHRGIN) .SIMTEMl IN,l) . SI MTEMl I N, 2 ) , S 1 MTEHl C N, 3 ) , TCOEFFI  N ) , 

3  THRSHLO(N) ,VEGTYPE(NI ) 

C  CHECK  THE  MANDATORY   ISOTHERMAL  DATE 

IFIYYMMDD  -    I SOTHRM t N , I YR ) )  20,10,20 

C  TURN  OFF   THE  DIFFUSION  MODEL,   SET  THE   PACK   TO  0  C   AND  ADJUST  THE 

C   ENERGY  BALANCE  ACCORDINGLY 

10  DREADYIN)   =  -1 

ENGBALINI   =  ENGBALIN)   ♦  CALDEF(N) 

RADLWN  =  RAOLWN  *  CALDEF(N) 

CALDEF(N)   =  0.0 

C  WEIGHT  AND  STORE  THIS  DATA 

20  WT  =  WEIGHT(N) 

COMPS(l)   =   COMPSll)    ♦   (TEMPMAX   »  WT ) 

C0MPS(2I   =   C0MPS(2)   ♦   ITEMPMIN  »  WT) 

C0MPSt3)   =  C0MPS13)   +   (PRECIP  •  WT) 

IF(PREWEOV(N)  1  '•0,'.0,30 
30  COMPSC)   =  COMPSU)   *   (ONTREESIN)  *  WT) 

C0MPSI5)   =  C0MPS(5)   ♦   (RADSWN  •  WT) 

C0MPS(6)   "  C0MPS(6)   *   (RADLWN  »  WT) 

C0KPSI7)   =  C0MPS(7)   *    lENGBAL(N)    *  WT) 

C0MPSI8)   =  C0MPS(8)   ♦   (CALDEFIN)   •  WT ) 

C0MPS(9)    =  C0MPS(9)   ♦    (PREWEQVIN)   »  WT) 
40  COMPSIIO)   =  COHPS(IO)   ♦   (WATERIN  •  WT) 

COMPS(ll)   =  COMPSdll   ♦   (EVAPOTR  •  WT) 

C0HPS(13)   =  C0MPS113)   ♦   (RECHRGIN)   »  WT) 

COMPSIm   =  COMPSIU)   ♦   IGENRO  •  WT ) 

CALL  ETCODE   ( ETFROH.COMPS ( 12 ) ) 

RETURN 

END 


Function  ROUNDF 

FUNCTION  ROUNDF  (SIGN) 
C  DETERMINE  THE  SIGN  FOR  ROUNDING 

IF(SIGN)  10,20,20 
10  ROUNDF  =  -0.5 

RETURN 
20  ROUNDF  =  0.5 

RETURN 

END 


Subroutine  WRITES 

SUBROUTINE  WRITE3  IK.COMPS) 

C  THIS  VERSION  OF  THE  OUTPUT  ROUTINE  PRINTS  ONE  COMPOSITE  LINE  PER 

C  DAY 

COMMON  AIRTEHC(25,6) 

COMMON  CALDEFI25) ,CDMAX 1 25 ) ,COVDEN(  25 ) 
COMMON  DREAOYI25) 

COMMON  ENGBAL(25) ,ET,ETDAILY(25,I2) 

COMMON  FREEWAT(25) 

COMMON  IS0THRMI25,20) 

COMMON  LASTUSDI25).LEVEL1,LEVEL2 

COMMON  MMDD 

COMMON  NDAYSN0I25) ,NDIVSBL,NSUB,NYEARS 
COMMON  0NTREESI25) 

COMMON  PEAKPPTI25,20) ,PEAKWEI25,20) ,PHASE (25) ,POTENT( 24) , 
I  PREWEQVI25) 
COMMON  RECHRG(25) 

COMMON  SIHTEMl(25,3),SUBIDI25,6),SLPASPI25,2'i) 

COMMON  TCOEFFI 25) ,THRSHLD(25) ,TOPLOT( ID 

COMMON  VEGTYPE (25) 

COMMON  WEIGHT(25) ,WSHEDID(6) 

COMMON  YEARSI20),YYMMDD 

INTEGER  DREADY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTYPE 

INTEGER  WSHEOID 

INTEGER  YEARS, YYMMOO 

COMMON/FOROATA/  FOOTNOT ( 26 ) .MAXLINE 
INTEGER  FOOTNOT 

C0MM0N/MASTER/0ATEI3) ,TMXMSTR,TMNMSTR,PPTMSTB,PPTONOW,OBSHYDR, 


I   POTRAD.MSTREOF , lYR 
INTEGER  DATE 

C OHMON / WATRB AL/ETF ROM, EVAPOTR,GENRO, PRECIP, RAO  IN, RAOLWN, RADSWN, 
1    TEMPMAX, TEMPMIN, WATERIN 
DIMENSION  COMPS I  16) ,FROM( 7) , I  TEMP  SI  3) ,OUT( 13) ,0UT1 ( 7) 
EOUIVALENCE    (OUT(  I ),OUTl I  1 ) ) 
DATA  L1NES/-1/ 

DATA  FR0M/3HC     ,3H  S   ,3HCS   ,3H     E,3HC  E, 3H  SE,3HCSE/ 

J   =  0 

IF(C0MP5(14).LE.0.0)   J  =  -1 
ITEHPS(l)   =  COMPSll)   ♦  ROUNDF  ICOHPSID) 
ITEMPSI2)   =  C0MPSI2)  ♦  ROUNDF  IC0MPS(2)) 
ITEHPS(3)   =  COHPSdS)   ♦  ROUNDF(C0MPS(15)  ) 
OUT( 1)    =  C0HPS(3) 
DUT(2)    =  C0MPS(16) 
IFICOMPSO)  )  10,10,50 
10  0UT1(3I   =  COHPS(IO) 
OUTICI    =  COMPSIU) 
OUTl (51    =   FR0M(C0MPS(12) ) 
0UT1(61    =  COMPS( 13) 
OUTl ( 7)    =  COMPSI 14) 
J   =   J   ♦  7 

C  DETERMINE  HOW   TO  WRITE   THE  LINE 

IF IK.EO. 1H$ ]   GO  TO  20 

WRITE    16.910)   K, I  TEMPS . (OUTl ( I ) ,1  =  1 ,J ) 
910  FORMAT ( 1XA10,2X3I4,2F6.2,50XF5.2,2XF6.4,2XA3,5XF5.2,6XF5.2) 

LINES   =   LINES   -  1 
RETURN 

C  CHECK  THE   LINE  COUNTER 

20    IF(LINES)  30,30,40 

30  WRITE    16,920)    FOOTNOT, WSHEOID, NSUB 
920  FORHATI«0«13A10/1X13AIO/»1*54X*WATER  BALANCE  SIMULATION*/ 

1  1X6A10,41X<'COMPOSITE   aF»I3,»  SUBSTATIONS'/ 

2  «0»l5x»TEMP   (F)       PRECIP   (IN)      INTERCEPTED     RAO  ICAL)   ENG  BAL 

3  SN0WP4CK  INPUT     EVAPOTRANS         RECHARGE  GENERATED*/ 

4  5X*DATE  MAX  MIN  AVE  DAY  ACCOM  (IN)  SW  LW  (CA 
5L)       TEMP(C)   WE(IN)        UN)      UN)     FROM         REO   (IN)     RUNOFF  (IN)*) 

LINES   =  MAXLINE 
40  WRITE   (6,930)   DATE , I  TEMP S, ( OUT  1 ( I ) , I  =  I , J ) 
930  FORMAT (*  *  3  I  3 , 3X3 1  4 , 2F6 . 2 , 50XF5 . 2 , 2 XF6 .4 , 2X A3 , 5XF5 . 2 , 6XF5 . 2 ) 
LINES   =   LINES  -  1 
RETURN 

C  ALL  THE   INFORMATION   IS  TO  BE  PRINTED 

50  00  60   I    =  3,13 

60  OUTl I )    =  COMPSI 1*1 ) 

OUT(ll)   =  FR0H(COMPSI 12) ) 
J  =  J   ♦  13 

C  DETERMINE  HOW  TO  WRITE  THE  LINE 

IF(K.E0.1H$)   GO  TO  70 

WRITE   (6,940)   K, I  TEMPS, (OUTl I ), 1  =  1, J) 
940  FORMAT ( 1 XA 1 0 , 2 X3 1  4 , 2F6 . 2 , 5 XF6.4 , 3X3F6 . 1 , 4X F4 . 1 , 2XF6 . 2 , 2XF5 . 2 , 2X 
I  F6.4,2XA3,5XF5.2,6XF5.2) 
LINES  =  LINES  -  1 
RETURN 

C  CHECK  THE  LINE  COUNTER 

70  IF(LINES)  80.80,90 

80  WRITE    16.920)   FOOTNOT . WSHED I  0, N SUB 

LINES  =  MAXLINE 
90  WRITE   (6.950)   DATE . I  TEMPS, lOUTI I ) , I  =  1,J) 
950  FORMA Tl*  *3I3,3X3I4,2F6.2,5XF6.4,3X3F6.l,4XF4.1,2XF6.2,2XF5.2,2X 
1  F6.4,2XA3,5XF5.2,6XF5.2) 
LINES   =   LINES  -  I 
RETURN 

C 

C  YEARLY  TOTALS 

C 

ENTRY  WRITOT 

C  DETERMINE  HOW  TO  WRITE  THE  LINE 

IFIK.EO.lHt)   GO  TO  100 
WRITE    (6,960)    K, (COMPSI I ), 1=1,5) 
96  0  FORMAT ( 1XA10,20XF6.2,49XF5.2,F8.4,9XF6.2,4XF7.2/ I 
LINES  =  LINES  -  2 
RETURN 

C  CHECK  THE  LINE  COUNTER 

100   IFILINES)  110,110,120 

110  WRITE   (6,920)   FOOTNOT, WSHEOID, NSUB 

LINES   =  MAXLINE 
120  WRITE   (6,9701   OA TE . ( COMPS ( I ) , I = I , 5 ) 

970  FORHAT(*OTOTALS  THR0UGH*3 I  3 , 7XF6. 2 , 49XF6. 2 , F8 . 4, *      ICHNG  =«F6.2,*) 
1*3XF7.2/ ) 
LINES  =  LINES  -  3 
RETURN 
END 

Program  DAILYO 

OVERLAY  (0LAYS,3,1] 
PROGRAM  DAILYO 

C  DAILY  COMPOSITE  OUTPUT,   NO  ALTERNATIVES 

COMMON  AIRTEMCI25,5) 

COMMON  CALDEFI25),CDMAXI25).COVDENI25) 
COMMON  DREA0YI25) 

COMMON  ENGBALI25) ,ET.ETDAILYI25.12I 

COMMON  FREEWAT(25) 

COMMON  !SOTHRM(25,20) 

COMMON  LASTUSD(25) .LEVEL1,LEVEL2 

COMMON  MMDD 

COMMON  NOAYSNOI 25) , NDI VSBL , NSUB , NYE ARS 
COMMON  ONTREESI25) 

COMMON  PEAKPPT(25,20) ,PEAKWE(25,20) ,PHASE(25).P0TENTI24), 
I  PREWEQVI25) 
COMMON  RECHRG(25) 

COMMON  SIMTEM1(25,3) , SUBI 0 1  25 , 5 ) , SLPAS P 1 25 ,  24 ) 
COMMON   TCOEFF (25) ,THRSHL0(25) ,T0PLOTI  11 ) 
COMMON  VEGTYPEI25) 
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COMMON  KEIGHK 25) .HSHEDID(5) 

COMMON   YEARS  I  20) . YYMHDD 

INTEGER  DRE4DY 

INTEGER  PHASE 

INTEGER  SUBID 

INTEGER  TOPLOT 

INTEGER  VEGTyPE 

INTEGER  HSHEOIO 

INTEGER  YEARS, YYMMOO 

CQMM0N/CMPSIT0/C0MPS(16) ,YRT0T(5) 
COMMON/ FORD AT A/   F OOT NOT  1 26  I , MAX L I NE 
INTEGER  FOOTnOT 

COMMON/MASTER/DATEI 3) .TMXMSTR,TMNMSTR,PPTMSTR,PPTON0W,0BSHYDR, 
1  POTRAD.MSTREOF.IYR 
INTEGER  DATE 
COMMON/ PLOTS/PLOT!  11  I 

COMMON /WaTRBAL/eTF ROM, EVAPOTR.GENRO.PRECIP, RAD  I N.RAOLHN.RAOSHN, 

1  temphax,tEhpmin,waterin 
integer  f00tntei26) 

data  f00tnte/260hn0rmal  simulation  only 

I 
2 
3 

C  DEFINE  THE  FOOTNOTE 

00  1  N  =  1,25 
1   FOOTNOTINI    =  FOOTNTE(N) 
C  MAXLINE  =  52  -  NUMBER  OF  ALTERATIONS 

HAXLINE  =  52 
C  WRITE  THE  FOOTNOTE  ON  THE  PLOT 

WRITE   (11,900)  FOOTNOT 
900  FORMAT! /•0»13A10/1X13A10//) 

C  READ  A  MASTER  CARD 

10  CALL  RDMSTR 

IF(MSTREDF)  20,20,130 

C  GENERATE   THE  DATA  AND  PERFORM  THE  SIMULATIONS  FOR  EACH  SUBSTATION 

20  DO  30  N  =  l.NSUB 

CALL   GENOATA  (N) 

CALL   N0RM3  (N) 
30  CONTINUE 
C  WRITE  OUT  THE  COMPOSITE  LINE 

COMPSdSl   =   (COMPSIl)   ♦  C0MPS(2))   »  0.5 

C0MPS(16)   =  C0MPS116)   ♦  C0HPS(3) 
C  CONVERT  THE  CALORIE  DEFICIT  TO  A  PACK  TEMPERATURE 

IFIC0HPSI9).NE.0.0)   COMPS(B)   =  -  COMPS ( 8 ) / ( COMPSI 9 )   •  1.271 

CALL   WRITE3    IIHJ, COMPS) 
C  ADO  THESE   VALUES   INTO   THE   YEARLY  TOTALS 

YRT0T(2)   =  YRT0T(2)   ♦  COMPSllO) 

YRT0T(3)   =  YRT0T!3)   t  COMPS(ll) 

YRT0T(5)   -   YRT0T(5)    ♦  COMPSd*) 
C  STORE   THE    INFORMATION  FOR   THE  PLOTS 

PLOTIl)   =  OBSHYDR 


PL0T(2I   =  C0MPS(9) 

PLOTI 3 ]   =  COMPSI 101 

PLOT!*)   =  COMPSIllI 

PL0TI5)   =  C0MPSI13) 

PL0TI6)   =  COMPSIID 
f-  PERFORM  THE  PLOTS  BETWEEN  APRIL   1  AND  SEPTEMBER  30 

IF(MMOD  -   401)  Il0,90,<i0 
'.0   IF(930  -  HHOD)  110,50,100 

C  ON  9/30,   PRINT  THE  YEARLY  TOTALS.   ZERO  THE  ACCUMULATED  PRECIP  AND 

C   RESET   THE  DIFFUSION  MODEL  SWITCHES 

50  YRTOT(l)   =  C0MPS(16) 

YRTOTCi)   =  COMPSI  131  -  YRTOTdl 

CALL  WRITOT  llHt.YRTOT) 

YRTOTU)   =  C0MPSI13) 

YRT0TI2)   =  0.0 

YRT0TI3I   =  0.0 

YRT0TI5)   =  0.0 

C0MPSI16)    =■  0.0 

00  50  N  =  1,NSUB 
50   DREADYIN)    =  0 

C  CHECK  THE  YEARS  TO  BE  SURE  ALL  DECKS  STILL  CORRESPOND  TO  THE 

Q  PARAMETER  DECK 

IFIDATEI3)   -  YEARS(IYR))  70,80,70 
70  WRITE   16,910)    DATE  I  3 ) , I YR , I YR , YE AR S II YR ) 
910  FQRMATItlTHE  DATA  DECKS  ARE  OUT  OF  PHASE  WITH  THE  WATERSHED  PARAME 
ITER  DECK.     THE  WATER  YEAR  JUST  ENDING  (19*12, «l    IS  DECK  NUMBER«I3, 
2  *,»/*     BUT  SPECIFIED  CONDITIONS  CARD  NUMBER»I3,*  FOR  EACH  SUBSTAT 
3I0N   IS  FOR  19«12) 
CALL  EXIT 
80   lYR   =   lYR   *  I 
GO  TO  100 

C  ON  4/1,   WRITE  THE  ORDINATE  LINE 

90  WRITE  (11,920) 
920  F0RMAT(9X124I1H. ) ) 

C  PERFORM  THE  PLOT,   THEN  ZERO  THE  COMPOSITE  LOCATIONS 

100  CALL  PLOTTER 
110  DO   120  N  I  1,14 
120  COHPSIN)   =  0.0 
GO  TO  10 

C  ALL  CARDS  HAVE  BEEN  READ.     IF  THE  DECK  ENDED  ON  9/30.  ALREADY 

C   CAUSING  THE  WATER  YEAR  TOTALS  TO  PRINT,   JUST  END  THE  JOB.     BUT  ON 

C   ANY  OTHER   DAY,   PRINT   OUT   THE  YEAR  TO  DATE   TOTALS  HERE 

130    IFIMMDD  -  930)  140,150,140 
140  YRTOTIll   =  C0HPS(15I 

YRT0TI4)   =  PL0TI5)   -  YRTOTUI 
CALL   WRITOT  (1H»,YRT0T) 

C  WRITE  THE  FOOTNOTE  ON  THE  LAST  PAGE 

150  WRITE    (5,930)  FOOTNOT 
930  F0RMAT(»0«13A10/1X13A10) 

C  RETURN  TO  THE  MAIN  UVERLAY  FOR  NORMAL  TERMINATION 

END 
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